# Limits to the development of feed-forward structures in large recurrent neuronal networks

^{1}Functional Neural Circuits Group, Faculty of Biology, Albert-Ludwig University of Freiburg, Germany^{2}Bernstein Center Freiburg, Albert-Ludwig University of Freiburg, Germany^{3}RIKEN Brain Science Institute, Wako, Japan^{4}Institute of Neuroscience and Medicine (INM-6), Computational and Systems Neuroscience, Research Center Jülich, Germany^{5}RIKEN Computational Science Research Program, Wako, Japan

Spike-timing dependent plasticity (STDP) has traditionally been of great interest to theoreticians, as it seems to provide an answer to the question of how the brain can develop functional structure in response to repeated stimuli. However, despite this high level of interest, convincing demonstrations of this capacity in large, initially random networks have not been forthcoming. Such demonstrations as there are typically rely on constraining the problem artificially. Techniques include employing additional pruning mechanisms or STDP rules that enhance symmetry breaking, simulating networks with low connectivity that magnify competition between synapses, or combinations of the above. In this paper, we first review modeling choices that carry particularly high risks of producing non-generalizable results in the context of STDP in recurrent networks. We then develop a theory for the development of feed-forward structure in random networks and conclude that an unstable fixed point in the dynamics prevents the stable propagation of structure in recurrent networks with weight-dependent STDP. We demonstrate that the key predictions of the theory hold in large-scale simulations. The theory provides insight into the reasons why such development does not take place in unconstrained systems and enables us to identify biologically motivated candidate adaptations to the balanced random network model that might enable it.

## 1 Introduction

For several decades now it has been commonly assumed that functional structures in the brain develop by strengthening synapses between neurons that fire in a correlated fashion (Hebb, 1949). In particular, feed-forward structures in which neural activity is propagated as a wave from one pool to the next (Abeles, 1991) would seem to be favored by a synaptic plasticity dynamics which strengthens causally correlated connections and weakens acausally correlated connections. This key property of spike-timing dependent plasticity (STDP) was postulated theoretically (Gerstner et al., 1993) before it was observed experimentally on the timescale of 10 ms (Markram and Sakmann, 1995; Markram et al., 1997; Bi and Poo, 1998), although strengthening of causal correlations on the timescale of 100 ms had already been found (Gustafsson et al., 1987).

Such feed-forward structures, also known as synfire chains or braids, have been postulated to underlie experimentally observed precise spike-timing, for example in mammalian cortex (Eckhorn et al., 1988; Gray and Singer, 1989; Abeles et al., 1993; Prut et al., 1998; Ikegaya et al., 2004; Pulvermüller and Shtyrov, 2009) or songbird HVC (Hahnloser et al., 2002; Kozhevnikov and Fee, 2007). The propagation of waves of activity along synfire chains has been shown to be stable under quite general conditions (Diesmann et al., 1999; Gewaltig et al., 2001). Apart from being a natural candidate for the representation of serial activities, such as the sequential activation of muscles to generate a movement, Bienenstock (1995) conjectured that all cortical computation is carried out by the action and interaction of such structures. It has since been shown that co-ordinated activity between synfire chains can realize the phenomenon of compositionality (Hayon et al., 2004).

Given the great degree of interest in feed-forward structures, the experimentally confirmed existence of a plasticity rule that is apparently predestined to generate them, and the growing availability of both tools suitable for large-scale neuronal network simulations (e.g., Morrison et al., 2005; Migliore et al., 2006; Gewaltig and Diesmann, 2007; Pecevski et al., 2009) and the high-performance computers on which to run them, it seems curious that there have been no truly convincing studies demonstrating their development.

It is in the nature of numerical studies that they use simplified models and therefore run the risk of generating artifacts, i.e., behavior that does not occur if a simplified model is replaced with a more realistic one. When investigating STDP in recurrent networks, there are at least two commonly chosen model simplifications that run a high risk of yielding such non-generalizable behavior. So far, all reports of structure development in a recurrent network have been based on network models which represent only a small fraction of the number of synaptic inputs a neuron typically receives (Hertz and Prügel-Bennet, 1996; Levy et al., 2001; Izhikevich et al., 2004; Iglesias et al., 2005; Doursat and Bienenstock, 2006; Jun and Jin, 2007; Masuda and Kori, 2007; Hosaka et al., 2008; Liu and Buonomano, 2009; Fiete et al., 2010). This is a serious issue, as STDP is driven by correlation between pre- and post-synaptic neurons. Scaling down the network size and increasing the strength of individual synapses increases the correlation between pre- and post-synaptic neurons. This leads to a strong competition effect between the inputs which can result in a very small number of inputs driving the post-synaptic neuron in a non-biological winner-takes-all manner. This is particularly so if, as for the majority of previous studies, an additive model for STDP is assumed, which has a strong symmetry breaking tendency leading to a bimodal distribution of synaptic strengths (Song et al., 2000). However, even very early experimental findings on STDP revealed that the strength of potentiation and depression is dependent on the initial strength of the synapse (Bi and Poo, 1998, but see also Debanne et al., 1996, 1999; Montgomery et al., 2001; Wang et al., 2005). Theoretical studies have revealed that even quite a small dependence on the weight leads to a qualitatively different behavior and a unimodal distribution of synaptic strengths (van Rossum et al., 2000; Rubin et al., 2001; Gütig et al., 2003; see also Morrison et al., 2008) which is more similar to the distributions observed experimentally (Turrigiano et al., 1998; Sjöström et al., 2001; Song et al., 2005). Moreover, weight-dependent “soft-bounded” rules can be more plausibly generated with a kinetic modeling approach than an additive “hard-bounded” rule (Zou and Destexhe, 2007). It is therefore not clear, and should not be assumed, that behavior observed in networks with low numbers of synapses per neuron and/or additive STDP are representative for networks with biologically realistic numbers of synapses and weight-dependent STDP. Indeed, a network exhibiting cortical connectivity levels and using a model for STDP with weight dependence fitted to experimental findings did not develop structure either spontaneously or as a result of repeated synchronous stimuli to a subset of neurons (Morrison et al., 2007).

One reason for this may have been that a critical additional mechanism is required. Many of the above mentioned studies include homeostatic regulating mechanisms such as normalization of weights on the axon and dendrite (Doursat and Bienenstock, 2006; Fiete et al., 2010), pre-synaptic activity dependent weight dynamics to maintain a given post-synaptic rate (Liu and Buonomano, 2009) or axonal pruning of weak synapses once a certain number of very strong synapses have been created (Jun and Jin, 2007). Additionally, STDP has been supplemented with plasticity dynamics on other time scales: Izhikevich et al. (2004) employ both short-term plasticity and slow long-term potentiation, whereas Jun and Jin (2007) implement slow long-term depression and a reversible activation and silencing of synapses. Clearly, some of these mechanisms are more biologically plausible than others and it may well turn out that one or more of them is necessary for the development of structure in networks with realistic connectivity.

In this article, we develop a theory for the recruitment into a structure of neurons in a recurrent network in response to repeated external synchronous input. We show that recruitment is characterized by an unstable fixed point that cannot be stabilized by simply introducing pre- or post-synaptic homeostatic mechanisms such as those proposed in previous studies (Doursat and Bienenstock, 2006; Fiete et al., 2010). We demonstrate that the key predictions of the theory can be reproduced by a large-scale neuronal network model. Finally, we consider what biologically plausible adaptations to the network model could allow the stable propagation of feed-forward structure.

## 2 Materials and Methods

### 2.1 Network Model

Our network model is based on the balanced random network model of Brunel (2000), see Figure 1A. The neurons are 80% excitatory and 20% inhibitory and randomly connected; excitatory–excitatory connections are subject to weight-dependent STDP, all other connections are static. The network activity is presumed to be in the asynchronous irregular regime with an average firing rate of *ν*_{r}.

**Figure 1. Network and single neuron model.** **(A)** Balanced recurrent network model based on Brunel (2000). A subset of the excitatory population (E_{stim}, pink disks) receives feed-forward input from each member of an external group of size *N*_{SIP} (black disks). The external group fires synchronously with rate *ν*_{s} and asynchronously with rate *ν*_{a}. The externally stimulated population E_{stim} makes recurrent connections to the rest of the excitatory population. *N*_{hc} neurons receive a number of connections from E_{stim} above a given threshold *C*_{hc} (high-connectivity group, purple disks) whereas all other neurons receive fewer connections from E_{stim} (blue disks). All neurons also receive external excitatory Poisson input (not shown). **(B)** Determination of the high-connectivity group from the connectivity matrix (shown here in diluted form). The subdivisions of the connectivity matrix indicate the source and target populations, e.g., *EI* contains the connections from the inhibitory population *I* to the excitatory population *E*. The entries in row *i* of the connectivity matrix indicate the sources of the synaptic inputs of neuron *i* according to the color scheme in **(A)**. The histogram gives the number of connections each neuron receives from E_{stim}. If neuron *i* receives a number of connections from E_{stim} greater than or equal to the arbitrary threshold *C*_{hc} (vertical dashed line) it is considered a member of the high-connectivity group *E*_{hc}. **(C)** Simplified input structure of a neuron in a recurrent network: the post-synaptic neuron *i* (gray disk) receives input from neurons from E_{stim} (pink disks) that fire asynchronously with rate *ν*_{a} and synchronously with rate *ν*_{s} and from neurons selected from the rest of the network (blue disks) that fire asynchronously with rate *ν*_{r}. The mean synaptic weights from these populations are *w*_{s} and *w*_{r}, respectively.

A subset E_{stim} of the excitatory population receives feed-forward input from each element of an external group of size *N*_{SIP}; these connections are also subject to STDP. Each neuron in the external group spikes independently as a Poisson process of rate *ν*_{a} and synchronously with the rest of the group as a Poisson process of rate *ν*_{s}, thus implementing a single interaction process (SIP, see Kuhn et al., 2003). The total rate of a neuron in the external group is the same as the network firing rate, i.e., *ν*_{a} + *ν*_{s} = *ν*_{r}. All neurons in the network receive additional excitatory independent Poisson spike trains; the neurons in the externally stimulated group receive Poisson spike trains at a slightly reduced rate to compensate for their additional external stimulus.

Due to the random connectivity, the number of connections neurons receive from E_{stim} is binomially distributed. For a given connectivity threshold *C*_{hc}, a high-connectivity group *E*_{hc} can be located within the excitatory population (not including E_{stim}). Each neuron in *E*_{hc} receives at least *C*_{hc} connections from E_{stim}, whereas all other excitatory neurons receive fewer than *C*_{hc} connections from E_{stim}. This is illustrated in Figure 1B. In general, *N*_{SIP}, *N*_{stim}, and *C*_{hc} can all be chosen independently. For a given binomial distribution of connections from E_{stim}, the choice of *N*_{stim} and *C*_{hc} completely determines *N*_{hc}, the number of neurons in the high-connectivity group *E*_{hc}. As we are investigating the development of feed-forward structures in this article, we are interested in the case that successive groups have similar sizes. We therefore link the independent variables such that *C*_{hc} = *N*_{SIP} and require *N*_{hc} ≈ *N*_{stim}. In other words, for a given value of *N*_{SIP}, we select the size of the externally stimulated population *N*_{stim} such that the binomial distribution of connections results in *N*_{hc} ≈ *N*_{stim} neurons that receive at least *N*_{SIP} connections from E_{stim} and all other excitatory neurons in the network receive fewer than *N*_{SIP} connections from E_{stim}.

### 2.2 Fixed Point Analysis of Structural Development

If pre- and post-synaptic spike trains are stochastic, the weight updates of a synapse can be described as a random walk. Using Fokker–Planck mean field theory, the drift of the random walk corresponds to the average rate of change of synaptic strength. The drift can be calculated as a function of the correlation of the pre- and post-synaptic spike trains. Assuming stationarity, the raw cross-correlation function of the pre-synaptic spike train and the post-synaptic spike train with mean firing rates of *ν _{i/j}* = 〈

*ρ*〉 is given by:

_{i/j}Γ_{ji} (Δ*t*) = 〈*ρ _{j}* (

*t*)

*ρ*(

_{i}*t*+ Δ

*t*)〉

_{t}

This results in the following expression for the mean synaptic drift:

where *F*_{±} (*w*) describe the weight-dependent potentiation/depression of a synapse due to a single spike pair and *K*_{±} (Δ*t*) = exp (−|Δ*t*|/*τ*_{±}) is the window function of STDP. A more thorough derivation of (1) can be found in Kempter et al. (1999) and Kistler and van Hemmen (2000; but see also Kempter et al., 2001; Rubin et al., 2001; Gütig et al., 2003; Morrison et al., 2008; Gilson et al., 2009a).

To investigate the fixed point behavior of STDP in the context of a large recurrent network in which structure development is being induced by synchronous stimulation of a subset of neurons, we make the following assumptions, abstracted from the network model displayed in Figure 1A. Each neuron *i* receives excitatory connections from neurons selected from the subset of neurons E_{stim} receiving synchronous stimulation and connections from neurons selected from the rest of the network. Synchronous spikes from the stimulated subset occur irregularly and are governed by a Poisson process with a rate of *ν*_{s}. Each neuron in the stimulated subset also fires independently as a Poisson process with rate *ν*_{a}; the composite activity for a neuron from the stimulated subset is therefore a Poisson process with rate *ν*_{a} + *ν*_{s}. The activity of a neuron in the rest of the network can be described as a Poisson process with rate *ν*_{r}. The input structure of an excitatory neuron in the recurrent network is illustrated in Figure 1C. We further assume that the post-synaptic spike rate is generated by an inhomogeneous Poisson process with scaling factor *γ*, threshold *ν*_{0} and membrane potential *u _{i}*, where [·]

_{+}denotes a piecewise linear function that is defined as 0 for negative values, i.e., [

*x*]

_{+}= (

*x*+ |

*x*|)/2. In the following, we will assume the function is positive and omit the square brackets. The membrane potential is given by where

*w*is the weight of the synapse between neuron

_{ij}*j*and neuron

*i*, and

*∊*(

*t*) describes the time course of the response of the membrane potential to the arrival of an excitatory event, i.e.,

*w*

_{ij}*∊*(

*t*) is the post-synaptic potential. Finally, we assume that the dynamics of the system is sufficiently well captured by considering the changes to the mean values of the synapses between the two input populations and the output neuron, and The expected rate of the output neuron is therefore:

where is the total area under the respective post-synaptic potential. The conditional firing rate after a spike at from a neuron in the unstimulated input population is:

Similarly, the conditional firing rate after a spike at from a neuron in the population receiving synchronous stimulus is:

Inserting these pre- and post-synaptic spike trains into the expression for the synaptic drift given in (1) results in a coupled system of differential equations for the evolution of the mean weights:

where The nullclines for *w*_{r} and *w*_{s} can be found by setting the left hand side of these equations to zero and solving for *w*_{r} and *w*_{s} respectively. The crossing points of the nullclines (if any) represent the fixed points of the system and can be found numerically.

To do this, appropriate values must be determined for the model parameters. For the numerical simulations in this study (see Section 2.3) we use a current-based integrate-and-fire neuron model with *α*-shaped post-synaptic currents, i.e., the dynamics of the membrane potential is given by:

where *τ*_{m} is the time constant of the membrane potential, *C* is the membrane capacitance and *I* is the input current. The time course of a post-synaptic current with peak amplitude *w* is given by:

We can therefore calculate the time course of the post-synaptic potential as:

and determine the model parameters.

and

in terms of the parameters of the integrate-and-fire neuron model and the STDP update rule.

### 2.3 Numerical Simulations

We perform numerical simulations of the network model described in Section 2.1. Our realization of the abstract model consists of 80,000 leaky integrate-and-fire neurons with *α*-shaped post-synaptic currents. The neurons are 80% excitatory and 20% inhibitory, each receiving 3840 excitatory and 960 inhibitory recurrent connections. Excitatory–excitatory connections are drawn from an initial normal distribution and are subject to weight-dependent power-law STDP (Morrison et al., 2007; Standage et al., 2007).

Each neuron except the externally stimulated group E_{stim} receives an additional independent Poisson spike train equivalent to 1200 excitatory inputs each spiking at *ν*_{x}. When no external synchronous stimulation is provided (*N*_{stim} = 0), a choice for the rate of the independent Poisson input to each neuron of *ν*_{x} = 14 Hz results in network activity in the asynchronous irregular regime (Brunel, 2000) with an average firing rate of *ν*_{r} and a unimodal distribution of the plastic recurrent weights with mean 38.6 pA and SD 3 pA (data not shown). In the absence of homeostasis, *ν*_{r} = 9.6 Hz; in the presence of weight normalizing pre-synaptic homeostasis applied every 3 s, *ν*_{r} = 9.5 Hz; in the presence of pre- and post-synaptic homeostasis applied alternately every 1.5 s, *ν*_{r} = 9.2 Hz.

When external synchronous stimulation is provided (*N*_{stim} > 0), each neuron in E_{stim} receives additional external stimulus at a total rate of *ν*_{r} from each of the *N*_{SIP} neurons in the external group E_{SIP}. Consequently, the total rate of the independent Poisson spike trains received by neurons in E_{stim} is reduced to 1200·*ν*_{x} − *N*_{SIP} *ν*_{r} to compensate for the additional external stimulus. We draw the initial weights of the recurrent plastic synapses in the stimulated network from a normal distribution with the same mean and SD as the equilibrium weight distribution of the unstimulated network. The initial weights of the connections from the external group E_{SIP} to the stimulated group E_{stim} are set to the mean of the equilibrium weight distribution of the unstimulated network.

During an initial equilibration period of 100 s the external group E_{SIP} generates independent Poisson spike trains at rate *ν*_{a} = *ν*_{r} with no synchronous events, i.e., *ν*_{s} = 0 Hz. After this period, the external group E_{SIP} provides partially synchronous stimulus to the stimulated group E_{stim} with a synchronous rate of *ν*_{s} = 3 Hz and an asynchronous rate of *ν*_{a} = *ν*_{r} − *ν*_{s}.

We perform simulations with different values for the size *N*_{SIP} of the external group. The threshold *C*_{hc} for membership of the high-connectivity group *E*_{hc} is set equal to *N*_{SIP}. For each value of *N*_{SIP}, we examine the connectivity of the network to determine empirically the value of *N*_{stim} that results in a high-connectivity group of approximately equal size, i.e., *N*_{hc} ≈ *N*_{stim}.

Following the standards suggested by Nordlie et al. (2009), the network model details are provided in tabular form in Table 1. Table 2 contains the values for the parameters used in the network simulations.

**Table 1. Tabular description of network model after Nordlie et al. (2009).**

The network simulations were carried out with a computation time step of 0.1 ms on the Stallo Linux cluster (Notur, UiT) and the JUGENE BlueGene/P supercomputer using the simulation software NEST (Gewaltig and Diesmann, 2007) at revision 8611. Simulation scripts are available from the authors on request.

### 2.4 Data Analysis

Spike times are recorded from the stimulated group, the high-connectivity group, and from 500 neurons randomly chosen from the excitatory population (random group). The times of the synchronous spikes of the external group are also recorded. Weights of all outgoing connections were recorded in intervals of 3 s from the external group, the stimulated group, the high-connectivity group and the random group.

The development of the mean weights shown in Figure 5 is determined from the outgoing weights of the external group and the weights of the connections from the stimulated group to the high-connectivity group. The weight distributions in Figures 6D,H are histograms (bin size 1 pA) of the weights of the connections between the specific neuronal groups recorded at 25 min.

To create the histograms in Figures 6A,C,E,G we record the number of spikes in each group in each 1 ms time bin within the 10 ms directly after a synchronous event over a period of 10 s and normalize by the number of neurons and the number of synchronous input events. Figures 6B,F and 7A–D show the development of the peri-stimulus spiking activity of the specific neuronal groups, which is the number of spikes within the 10 ms before (control) and within the 10 ms after each synchronous input event. Data is normalized with respect to the number of neurons in a group and the number of synchronous input events within successive periods of 10 s.

The scatter plots in Figure 7E are generated by recording the time-to-first-spike following each synchronous event in a 1 min period for all excitatory neurons in the network. Each marker in the plot shows the mean and SD ofthe response times for a specific neuron; the color of the markers indicates the likelihood of the neuron to respond to a synchronous stimulus within 10 ms.

## 3 Results

### 3.1 Fixed Point Analysis of Weight Dynamics

To investigate the dynamics of STDP in a recurrent network undergoing synchronous stimulation, we select the update rule proposed by van Rossum et al. (2000), which is additive for potentiation (*F*_{+}(*w*) = *λ*) and multiplicative for depression (*F*_{−}(*w*) = *αλw*). This rule is a good fit to the experimental data on the dependence of the magnitude of synaptic strength change on the initial synaptic strength (Bi and Poo, 1998). Furthermore, this choice has the advantage that the rule exhibits the same qualitative behavior as other weight-dependent rules (see Morrison et al., 2008), yet is more tractable than many other weight-dependent rules such as the partially multiplicative model proposed by Gütig et al. (2003) or the power-law formulation (Morrison et al., 2007; Standage et al., 2007).

Inserting *F*_{±}(*w*) into the equation system for the synaptic drifts given in (5) and (6) and setting the left hand sides to zero results in the following quadratic expressions for the nullclines of *w*_{r} and *w*_{s}:

The nullclines can be determined analytically; the crossing of the nullclines can be found numerically. To this end we assume input to the neuron analogous to that received by a neuron in the recurrent network simulation described in Section 2.3. The rate of the asynchronous background population is *ν*_{r} = 9.6 Hz, and the asynchronous and synchronous rates of the SIP input are *ν*_{a} = 6.6 Hz and *ν*_{s} = 3 Hz. The total number of incoming recurrent excitatory synapses is We set STDP parameters *α* = 0.0132 pA^{−1}, *λ* = 0.03 pA, *τ*_{+} = 15 ms and *τ*_{−} = 30 ms and neuronal parameters *γ* = 3.66 × 10^{3}(Vs)^{−1} and *ν*_{0} = 176.1 Hz. In the absence of synchronous stimulation this choice of synaptic and neuronal parameters results in a self-consistent output rate *ν _{i}* =

*ν*

_{r}and a stable fixed point thus reproducing the corresponding values of the unstimulated network simulation. The parameters and are calculated from the neuronal and STDP parameters as described in Section 2.2. The nullclines for are illustrated in Figure 2A. Note that as the nullcline for

*w*

_{s}is much more sensitive to the value of

*w*

_{r}than vice versa. The value of is negative for values of

*w*

_{s}above the

*w*

_{s}-nullcline and positive for values of

*w*

_{s}below it. Similarly, is positive for values of

*w*

_{r}to the left of the

*w*

_{r}-nullcline and negative for values to the right of it. Consequently, the fixed point indicated by the intersection of the nullclines is a stable attractor. Figure 2B shows the location of the fixed point as a function of the number of synchronous inputs . Again,

*w*

_{s}is the more sensitive variable.

**Figure 2. Fixed point of the synaptic drift equations.** **(A)** Nullcline for *w*_{s} as a function of *w*_{r} (red curve) and for *w*_{r} as a function of *w*_{s} (blue curve) for . The crossing point of the nullclines determines the fixed point of the coupled equations. Arrows show the sign of the derivatives of *w*_{r} and *w*_{s}; arrows pointing up (down) indicate arrows pointing right (left) indicate for the values of *w*_{r} and *w*_{s} at the base of the arrows. **(B)** Location of fixed point (; blue, ; red) as a function of the number of synchronous inputs .

We can now exploit our knowledge of the network connectivity to determine the distribution of and thus the distribution of fixed points. Each neuron receives *C*_{E} connections drawn randomly from the *N*_{E} excitatory neurons in the network, so the number of connections drawn from the E_{stim} population receiving synchronous stimulus is distributed binomially; the number of trials is *C*_{E} and the probability of success in an individual trial is *p* =*N*_{stim}/*N*_{E}, where *N*_{stim} is the number of neurons in E_{stim}. Multiplying the probability of drawing *C*_{s} connections P_{P} (*C*_{s}, *C*_{E}) with the number of excitatory neurons *N*_{E} gives the expected number of neurons in the network that have *C*_{s} connections from *E*_{stim}. Let us now consider an idealized network model in which the number of neurons with *C*_{s} connections from E_{stim} is indeed ⌊*N*_{E} · *P _{p}* (

*C*

_{s},

*C*

_{E})⌋ and assign each neuron a connectivity rank

*R*, where

^{i}*R*<

^{j}*R*entails Figure 3A shows the equilibrium average strength of the synapses from the stimulated neurons as a function of the size of the stimulated group

^{k}*N*

_{stim}and the connectivity rank of the post-synaptic neuron in the idealized network. Only the neurons with the highest connectivity are displayed (

*R*≤ 1000), as these are the neurons that will most easily be recruited into a structure. However, Figure 3A suggests that this is unlikely to occur, as increases slowly with decreasing

^{i}*R*, i.e., neurons which receive similar numbers of connections from the stimulated group will tend to develop similar synaptic strengths. Consequently, there will be no development of clearly defined groups that respond in a qualitatively different fashion to a synchronous event in the stimulated group. The increase of with increasing

*N*

_{stim}is more rapid; small variations in the size of a synchronously firing group can have a large effect on the outgoing weight distribution.

**Figure 3. Distribution of fixed points and likelihood of responding to synchronous stimulus.** **(A)** Average strength of synapse from the stimulated group as a function of the size of the stimulated group (horizontal axis) and the connectivity rank of the post-synaptic neuron (vertical axis). The black line indicates the group of highest connectivity equal in size to the stimulated group (*R* = *N*_{stim}). **(B)** Expected number of spikes emitted by the post-synaptic neuron in the 15 ms following a synchronous spike in the stimulated group; axes and black line as in **(A)**; the black curve indicates the number of neurons which respond to a synchronous event with a spike with probability greater than or equal to the recruitment threshold Θ_{r} = 0.5. The black line and curve represent a return map for the development of feed-forward structure in a network with an unstable fixed point at *N*_{stim} = 701. The blue arrows show two iterations of the return map for an initial stimulated population of size *N*_{stim} = 697; the pink arrows show one iteration for *N*_{stim} = 706. **(C)** Location and nature of the fixed point as functions of the recruitment threshold Θ_{r}. Dashed curve: fixed point , solid curve: slope of the tangent to the *N*_{r} (*N*_{stim}) curve at .

This insight can be further clarified by calculating the number of spikes a post-synaptic neuron is expected to produce within a given time *t* after a synchronous event in the stimulated input group. This can be found by integrating the conditional firing rate:

Figure 3B shows the number of spikes expected within 15 ms after a synchronous event for every fixed point calculated in Figure 3A. To determine the development of functional structure, we introduce a recruitment threshold Θ_{r}. All neurons that have a probability of responding to a synchronous event in E_{stim }of at least Θ_{r} are considered to have been recruited into the structure, neurons for which *P*_{15} (spike) < Θ_{r} are considered to be outside the structure. The black curve indicates the number of neurons recruited *N*_{r} as a function of the size of the stimulated group *N*_{stim} for an arbitrary choice of recruitment threshold Θ_{r} = 0.5, i.e., for a choice of *N*_{stim}, the height of the black curve gives the number of neurons that are recruited by this stimulus such that each neuron has at least a 50% probability of firing immediately after a synchronous event in E_{stim}.

The black diagonal line indicates a connectivity rank equal to the size of the stimulated group, i.e., the *N*_{stim} neurons in the network that receive the most input from the synchronously stimulated group. The two lines represent an approximation of a return map for groups of synchronously firing neurons: the curve indicates *N*_{r} (*N*_{stim}) and the diagonal indicates *N*_{stim}, therefore their intersection reveals the fixed point at which the number of neurons recruited is equal to the number of neurons stimulated, If the fixed point were stable, the recruited neurons would recruit another group of neurons, which would in turn recruit another such group, thus stably propagating structure through the network. However, the fixed point is unstable, since the slope of the *N*_{r} (*N*_{stim}) curve is positive at the intersection. If *N*_{stim} is initially below the crossing point, the number of neurons recruited is smaller than *N*_{stim}. These neurons will in turn recruit even fewer neurons, and so the structure will decay to nothing within a few iterations. Conversely, if the number of stimulated neurons is larger than the crossing point, the number of neurons recruited will be greater than *N*_{stim}. These neurons will in turn recruit an even greater number of neurons, until the entire network is recruited within a few iterations. This is illustrated by the map iterations shown as arrows in Figure 3B. Note that the iterations are approximative, as the composition of pre- and post-synaptic rates is in general not constant over iterations.

Figure 3C shows the dependence of the fixed point on the choice of the recruitment threshold Θ_{r}. As the requirement on the reliability of spiking following a synchronous event increases, the fixed point shifts to greater values of *N*_{stim}. The slope of the tangent of the *N*_{r} (*N*_{stim}) curve is positive across the whole range of Θ_{r}, demonstrating that the fixed point is always unstable.

The main factors influencing the shape of Figure 3 are the random wiring of the network and the smooth change of the fixed point of the weight distribution with increasing number of synchronous connections (Figure 2B). The number of connections a neuron receives from a given group is binomially distributed. This causes a “fan-out” structure – if we can find a group of neurons which receives at least *k* inputs from a synchronously active group, we can also find an even larger group which receives slightly fewer than *k* connections. Due to the smooth change of the fixed point of the weight distribution with *k*, the equilibrium weight distribution and thus the equilibrium response to a synchronous event of a neuron receiving *k *−* *1 or *k *−* *2 synchronous connections is very similar to the equilibrium distribution and response of a neuron receiving *k* connections.

It is therefore worthwhile to determine to what extent the smooth shift of the weight distribution shown in Figure 2B is robust to parameter changes. Figure 4A shows the development of fixed points as a function of the number of synchronous connections for four different scenarios. In the top left panel we consider the case of lower input rates: *ν*_{r} = 5 Hz, *ν*_{a} = 3.5 Hz and *ν*_{s} = 1.5 Hz, corresponding to a lower recurrent network rate and a lower synchronous stimulus rate. We set *ν*_{0} = 190 Hz and *γ* 7.02 × 10^{3 }(Vs)^{−1}, which corresponds to a reduced external Poisson rate and results in a self-consistent rate *ν _{i}* =

*ν*

_{r}when there is no synchronous stimulus . The top right panel reproduces the standard input scenario (

*ν*

_{r}= 9.6 Hz,

*ν*

_{0}= 175 Hz) except the rate of synchronous events is a greater proportion of the total rate of the synchronous population:

*ν*

_{a}= 2.6 Hz and

*ν*

_{s}= 7 Hz. The bottom left panel reproduces the standard input scenario (

*ν*

_{r}= 9.6 Hz,

*ν*

_{a}= 6.6 Hz,

*ν*

_{s}= 3 Hz,

*ν*

_{0}= 175 Hz) for the multiplicative STDP model introduced by Rubin et al. (2001):

*F*

_{+}(

*w*) =

*λ*(

*w*

_{max}−

*w*) and F

_{−}(

*w*) =

*αλw*. The synaptic parameters are chosen as for the additive/multiplicative model of van Rossum et al. (2000) investigated above, but with

*α*= 2.13 and

*w*

_{max}= 200 pA. Finally, the bottom right panel shows the results for the standard input scenario in combination with the power-law STDP developed in Morrison et al. (2007) and investigated numerically in Section 3.2: and F

_{−}(

*w*) =

*λα*w. Synaptic parameters are set as above, but with

*α*= 0.057,

*w*

_{0}= 1 pA and

*μ*= 0.4. In this example, the nullclines cannot be calculated analytically. For both the alternative STDP models investigated, the deviating parameters are chosen to give the same results for the standard input scenario with no synchronous input as the additive/multiplicative model of van Rossum et al. (2000) investigated above: a self-consistent rate and a stable fixed point .

**Figure 4. Robustness of results with respect to input rate, STDP model and post-synaptic homeostasis.** **(A)** Location of fixed point (; blue, ; red) as a function of the number of synchronous inputs for different scenarios, see text for details. Top left: low recurrent and external rates (solid curves) with additional homeostatic regulation of the post-synaptic rate (dashed curves). Top right: high synchronous rate (solid curves) with additional homeostatic regulation of the incoming synaptic weights of the post-synaptic neuron. Bottom left: standard input rates and multiplicative STDP model. Bottom right: standard input rates and power-law STDP model. **(B)** Expected number of spikes emitted by the post-synaptic neuron in the 15 ms following a synchronous spike in the stimulated group in the low rate scenario; axes and black curves as in Figure 3. **(C)** As in **(B)** for the standard input rate scenario with the power-law STDP model.

It has previously been suggested that homeostatic regulatory mechanisms are crucial for the development of structure (Doursat and Bienenstock, 2006; Liu and Buonomano, 2009; Fiete et al., 2010). The equation system for the synaptic drifts given in (5) and (6) can easily be extended to incorporate post-synaptic homeostatic mechanisms. A term −*λ _{ν}* (

*ν*

_{i}−

*ν*

_{h}) acts to increase the synaptic weights if the rate of the post-synaptic neuron is less than a desired rate

*ν*

_{h}and to decrease them if the post-synaptic rate exceeds

*ν*

_{h}. A natural choice for

*ν*

_{h}is the self-consistent rate of the unstimulated network,

*ν*

_{r}. Similarly, a homeostatic mechanism to conserve the sum of a neuron’s incoming synaptic weights can be modeled by an additional term where is the total number of incoming synapses and

*w*

_{h}is a desired mean synaptic weight. A natural choice for

*w*

_{h}is the fixed point of the one-dimensional system in the case that the neuron receives no synchronous input . For both models of post-synaptic homeostasis, the expressions for the nullclines of

*w*

_{r}and

*w*

_{s}given in (12) and (13) acquire additional contributions to the linear and constant coefficients but maintain their quadratic form. In Figure 4A, the development of the fixed point assuming post-synaptic rate homeostasis is shown for the low rate scenario with

*λ*= 3 × 10

_{ν}^{−4 }pA. Similarly, the fixed point development assuming post-synaptic weight homeostasis is shown for the high synchronous rate scenario with

*λ*= 5 × 10

_{w}^{−6}s

^{−1}. In both cases, the presence of post-synaptic homeostasis causes a greater separation to develop between the mean weights of the unsynchronized inputs and the synchronized inputs.

In all scenarios, a smooth shift of the fixed point of the weight distribution with increasing number of synchronous connections can be observed. When combined with the assumption of random connectivity, all scenarios result in a qualitatively similar return map to that shown in Figure 3 with an unstable fixed point. This is demonstrated in Figures 4B,C with a recruitment threshold of Θ_{r} = 0.5 for the low rate scenario (without homeostasis) and the power-law STDP scenario, respectively. As the theoretical model is developed from the point of view of the post-synaptic neuron, it is not easily extensible to account for pre-synaptic homeostatic mechanisms such as the normalization of outgoing weights. However, due to the similarity in behavior between neurons with similar numbers of synchronous inputs, implementing a pre-synaptic mechanism to regulate the total strength of outgoing synapses will not result in symmetry breaking for some specific value of . These results lead us to conclude that stable development of structure as a response to synchronous stimulation cannot occur in randomly connected networks incorporating weight-dependent STDP without additional network, cellular or sub-cellular assumptions.

### 3.2 Simulation Results

To check the predictions of the theoretical model we simulate a large-scale neuronal network with biologically realistic numbers of incoming synapses per neuron and weight-dependent STDP, see Section 2.3 for details. An excitatory sub-population of the network E_{stim} of size *N*_{stim} is stimulated by an external group E_{SIP} of size *N*_{SIP} in which the neurons fire mostly independently but with synchronous events at 3 Hz. Figure 5 shows the development of the mean synaptic weight between the external synchronous group and the stimulated group (top panel) and between the stimulated group and the high-connectivity group (bottom panel) for different sizes of the external group *N*_{SIP}. The size of the stimulated group *N*_{stim} is chosen such that a high-connectivity group E_{hc} can be found of similar size (*N*_{hc} ≈ *N*_{stim}) where each neuron in the high-connectivity group receives at least *N*_{SIP} connections from the stimulated group whereas all other excitatory neurons in the network receive fewer than *N*_{SIP} connections from the stimulated group. For low values of *N*_{SIP}, the mean synaptic weight between E_{SIP} and E_{stim} saturates within 10 min. As predicted by the theoretical model derived above, the downstream effect is much weaker; little change can be seen in the mean synaptic weight between E_{stim} and E_{hc}, even though each neuron in the E_{hc} receives at least as many inputs from E_{stim} as each neuron in E_{stim} receives from E_{SIP}. For higher values of *N*_{SIP}, the network is still evolving at 10 min, but the substantially greater effect on the mean synaptic weights between E_{SIP} and E_{stim}, and between E_{stim} and E_{hc} can be clearly seen.

**Figure 5. Evolution of the mean synaptic weights as a function of the size of the external group N_{SIP}.** Light to dark curves correspond to

*N*

_{SIP}= 10, 20, 30 resulting in a size for the stimulated group

*N*

_{stim}= 46, 164, 295 respectively.

**(A)**Mean weight of synapses between the external group E

_{SIP}and the stimulated group E

_{stim}as a function of time.

**(B)**Mean weight of synapses between the stimulated group and the high-connectivity group receiving at least

*N*

_{SIP}connections from it (

*N*

_{hc}= 44, 158, 291).

The rapid change in network behavior as the number of synchronous inputs increases can be seen in Figure 6, see Section 2.4 for the data analysis details. In the upper panel, *N*_{SIP} = 20 and *N*_{stim} = 164. Each neuron in the stimulated group E_{stim} receives input from each of the 20 neurons in the external stimulating group E_{SIP}. Each neuron in the corresponding high-connectivity group E_{hc} of a similar size (*N*_{hc} = 158) receives at least 20 synaptic inputs from E_{stim}. All other excitatory neurons in the network receive 19 or fewer synaptic inputs from E_{stim} and on average 10. The probability that neurons in the stimulated group fire in response (i.e., within 10 ms) to a synchronous stimulation event increases slightly in the first 10 min and then saturates. The probability that neurons in E_{hc} fire in response to a synchronous stimulation event mediated by E_{stim} remains at base level (i.e., determined by the network firing rate of ≈10 Hz) throughout the simulation. The weight distributions between all groups are very similar (see Figure 6D).

**Figure 6. Evolution of network response to synchronous stimulus.** **(A)** Probability of emitting a spike after a synchronous event in the first 10 s of stimulation for the stimulated group (top), the high-connectivity group (middle) and the random group (bottom) for *N*_{SIP} = 20 and *N*_{stim} = 164. **(B)** Evolution of the probability of emitting a spike in the 10 ms immediately before (dashed lines) or immediately after a synchronous event (solid lines) in successive 10 s periods. Colors as in **(A)**. **(C)** As in **(A)** but for the 10 s period following 25 min of simulation. **(D)** Distribution of incoming synaptic weights at 25 min for the stimulated group (top) from the random group (blue), the synaptic weight distribution from the external group (mean 56.9 pA) is not shown; for the high-connectivity group (middle) from the random group (blue) and the stimulated group (pink); for the random group (bottom) from the random group (blue) and the stimulated group (pink). **(E–H)**: As for **(A–D)** but for *N*_{SIP} = 25 and *N*_{stim} = 228. The synaptic weight distribution from the external group to the stimulated group (mean 161.4 pA) is not shown in the top panel of **(H)**.

In the lower panel, *N*_{SIP} = 25, *N*_{stim} = 228 and *N*_{hc} = 228. Each neuron in E_{hc} receives 25 or more synaptic inputs from E_{stim}, whereas all other excitatory neurons in the network receive 24 or fewer inputs from E_{stim} and on average 14. The probability that neurons in E_{stim} fire in response to a synchronous stimulation event increases steadily over the course of the simulation until almost all neurons are firing in response to every synchronous event. The probability that neurons in E_{hc} fire in response to a synchronous stimulation event, mediated by E_{stim}, remains at approximately base level for the first 20 min. Between 20 and 25 min, the probability of response increases sharply. However, the probability that neurons randomly selected from the rest of the network fire in response to a synchronous stimulation event also increases sharply during this period. This is reflected in the weight distributions. The middle panel of Figure 6H shows the distribution of incoming synaptic weights to the high-connectivity group E_{hc}. The mean of the distribution of weights between E_{stim} and E_{hc} shifts to a higher value than that between randomly selected neurons and E_{hc}, as seen in Figure 2B. However, the same is true to a lesser extent for the distributions of the incoming weights to the random group (bottom panel of Figure 6H). Note that the mean of the distribution of the incoming synaptic weights to the stimulated group from the rest of the network has shifted to a lower value (top panel of Figure 6H) due to the acausal correlation between activity in the network and activity in the stimulated group. If the synapses from E_{SIP} to E_{stim} are clipped to an upper bound, this drift eventually leads to a decoupling of the stimulated group, as was reported in Morrison et al. (2007). However, the insight that STDP reduces the strength of recurrent connections was originally described by Abbott and Nelson (2000), and more recently proved for polysynaptic loops of arbitrary length by Kozloski and Cecchi (2010).

Between 27 and 28 min the increasing response of the network to synchronous events causes it to enter a pathological state in which some neurons are firing very rapidly all the time and some are almost silent. The theoretical model developed in Section 3.1 does not predict the occurrence or the characteristics of the pathological state, as it does not take into account the interactions between neurons in the network, only the interactions between a synchronously firing group and its post-synaptic targets.

To test the conclusion that the introduction of homeostatic mechanisms do not alter the nature of the fixed point, we repeat the above experiments with pre-synaptic homeostasis. Every 3 s the outgoing synapses of every excitatory neuron in the recurrent network are adjusted by a multiplicative factor such that the total sum of its plastic synapses is reset. Similarly to the results obtained without homeostasis, synchronous stimulation from a group of *N*_{SIP} = 20 external neurons is neither enough to produce a reliable response of the stimulated group (*N*_{stim} = 164) nor any substantial changes in the high-connectivity group (*N*_{hc} = 158), whereas a slightly stronger stimulus (*N*_{SIP} = 25, *N*_{stim} = 228, *N*_{hc} = 228) causes a network to enter a pathological state (see Figures 7A,B). The transition occurs much later, at around 90 min rather than 28 min in the absence of homeostasis, demonstrating the stabilizing effect of the weight normalization. Unlike the results shown in Figure 6F for pure STDP, as soon as the stimulated group reflects the stimulus with sufficient reliability, a clear difference between the response likelihood of the high-connectivity group and the network as a whole can be observed. The high response likelihood of the high-connectivity group breaks down rapidly when the network starts responding more strongly to the synchronous stimulation at around 60 min. This is another example of the decoupling of a synchronously active group from its embedding network as described above and in Morrison et al. (2007). As the response to the synchronous stimulus of the network generally lags behind the response in the high-connectivity group, acausal correlations are generated in the incoming synapses of the high-connectivity group. As the synaptic weights decrease, so does the mean membrane potential, thus diminishing the response of the neurons to the synchronous stimulus events. The stimulated group E_{stim} does not decouple in this fashion because the increasing strength of the connections from the external input can compensate for the loss of input from the recurrent network.

**Figure 7. Spike-timing dependent plasticity and homeostasis.** **(A,B)** Evolution of network response to synchronous stimulus where STDP synapses in the recurrent network are subject to weight normalizing pre-synaptic homeostasis for *N*_{SIP} = 20, *N*_{stim} = 164, *N*_{hc} = 158 and *N*_{SIP} = 25, *N*_{stim} = 228, *N*_{hc} = 228, respectively. Probability of emitting a spike in the 10 ms immediately before (dashed lines) or immediately after a synchronous event (solid lines) in successive 10 s periods for the stimulated group (pink), the high-connectivity group (purple) and the random group (blue). **(C,D)** As in **(A,B)** but for pre- and post-synaptic homeostasis; *N*_{SIP} = 50, *N*_{stim} = 584, *N*_{hc} = 576 and *N*_{SIP} = 75, *N*_{stim} = 961, *N*_{hc} = 963, respectively. **(E)** Response latency and reliability for a network with STDP, pre- and post-synaptic homeostasis. Mean and SD of the response interval for each excitatory neuron to each synchronous event in the 26th minute (left, middle) or 21st minute (right). Marker color indicates the likelihood of the neuron to fire within 10 ms of a synchronous event. From left to right the panels correspond to *N*_{SIP} = 20, 50, 75 and *N*_{stim} = 164, 584, 961.

Figures 7C,D show the development of the response likelihood when post-synaptic homeostasis is also applied, such that the total strength of incoming and outgoing synapses are normalized alternately every 1.5 s. The network behavior is much more stable. The network remains in a stable asynchronous irregular activity regime even when an external group of size *N*_{SIP} = 50 stimulates a group of *N*_{stim} = 584 excitatory neurons (Figure 7C). However, the response probability of the high-connectivity group (*N*_{hc} = 576) saturates at a lower level than the maximum response likelihood reached in Figure 7B with purely pre-synaptic homeostasis. Although increasing *N*_{SIP} increases the number of neurons in E_{stim} available to drive a potential high-connectivity group *E*_{hc}, it also increases the size of *E*_{hc}. In the presence of pre-synaptic homeostasis, this reduces the mean strength of the connections between *E*_{stim} and *E*_{hc}.

Synchronous stimulation from an external group of size *N*_{SIP} = 75 drives the network into a pathological state despite the homeostasis as can be seen in Figure 7D. The interaction between the homeostasis and the strong stimulus sets up an oscillatory behavior in which the network repeatedly enters a pathological state but is reset by the application of homeostasis. Between pathological interludes, all neurons have a high probability of responding to a synchronous stimulus. The absence of clearly defined neuronal groups can be seen in Figure 7E. In the left panel, the network with STDP, pre- and post-synaptic homeostasis shows little response to synchronous stimulation by a small group of *N*_{SIP} = 20 external neurons: the average latency is approximately the inter-spike interval of the network background rate. The middle and the left panel correspond to Figures 7C,D, i.e., synchronous stimulation by much larger external groups of *N*_{SIP} = 50 and *N*_{SIP} = 75 neurons, respectively. The neurons of the stimulated group (*N*_{stim} = 584 and *N*_{stim} = 961) can be seen to have separated themselves and respond quickly and reliably to the synchronous stimulus. The mean response latency shifts to an earlier value and the mean SD of the latency is reduced, demonstrating that the network as a whole responds earlier and more reliably to the stronger stimulus. Those neurons that respond earliest also have the lowest SD and are thus most likely to respond within 10 ms of a synchronous stimulus. However, there is no separation of a second group from the mass. We therefore conclude that no development of feed-forward structure can be observed, even in the presence of pre- and post-synaptic homeostasis and strong synchronous stimulus.

## 4 Discussion

We have developed a simple model using a mean field approach and a linear neuron model to investigate the propagation of feed-forward structure in plastic recurrent networks. The key prediction of the model is that the number of neurons recruited by a repeated synchronous stimulus protocol is subject to an unstable fixed point. A synchronously firing group of neurons of a size below that of the fixed point recruits a smaller group, leading to a failure of the structure to propagate, whereas a synchronously firing group of a size above that of the fixed point recruits a larger group, causing the whole network to be recruited. In other words, a repeated synchronous stimulus is always either not enough or too much. We demonstrated by simulation that a large-scale network behaves as predicted by the mean field theory.

The existence of the unstable fixed point is robust. Employing a different weight-dependent STDP model, altering the input rates or introducing homeostatic mechanisms, e.g., regulating incoming synapses so as to maintain the total synaptic strength or post-synaptic rate, only influences the dynamics quantitatively. The unstable fixed point can be shifted to a larger or smaller group size, but it cannot be turned into a stable fixed point. The instability arises through the combination of a binomial distribution of connections from a given synchronously active group with the smooth change of the stable fixed point in the synaptic weight distribution as a function of the number of synchronous connections. Thus, any change to the model assumptions which does not produce a qualitative change in at least one of these factors will produce similar results, e.g., distributing the conduction delays or defining the conduction delay to be largely axonal, although such adaptations of the network model may well suppress or alter the pathological state the network enters when the amount of stimulation is too high (see Lubenov and Siapas, 2008). As an example, we demonstrated that the introduction of pre- and post-synaptic homeostatic mechanisms that maintain the total strength of outgoing and incoming recurrent synapses makes the network more stable, but does not enable the development of feed-forward structures.

Stable propagation of feed-forward structure is only possible if one or more of our assumptions are false. Other simulation studies have previously reported structure development (Hertz and Prügel-Bennet, 1996; Levy et al., 2001; Izhikevich et al., 2004; Doursat and Bienenstock, 2006; Jun and Jin, 2007; Masuda and Kori, 2007; Hosaka et al., 2008; Liu and Buonomano, 2009; Fiete et al., 2010). However, it is our contention that the risk of non-generalizable results is high when studying networks that are small, have few incoming synapses per neuron, or assume STDP rules that are independent of the synaptic weight. Each neuron in our network model receives a biologically realistic number of inputs (6000) and the network exhibits a biologically realistic degree of sparseness (connection probability of 0.06), while our model of STDP reproduces the experimentally found weight dependence of the amount of potentiation and depression a synapse undergoes (Bi and Poo, 1998; Morrison et al., 2007). We therefore suggest that the previously reported development of structure be considered artifactual until a greater degree of validation can be obtained. In our view, a weak proof of principle would be the demonstration of bounded structural growth in a network model with at least as great a degree of biological realism as in our study with respect to connectivity and plasticity, with additional assumptions that at least do not contradict experimental findings. A stronger, more convincing, proof would be provided by such a network model in which the additional assumptions are experimentally motivated.

Given that the brain does develop stable signal pathways, there is every reason to expect that a network model can be developed that fulfills these criteria for a strong proof of principle. One indication that the model we have investigated is overly simplified is that the combination of network structure, neuronal dynamics and plasticity leads to a symmetrical distribution of synaptic weights, rather than the skewed distribution observed experimentally (Sjöström et al., 2001). Therefore, a promising candidate adaptation that might result in a stable fixed point and thus stable propagation of feed-forward structure is the initial connectivity of the network. Here, we assume a random graph; however, the cortical network consists of several layers with layer-specific connectivity and exhibits long range patchy connections (Lewis et al., 2002; Schüz and Braitenberg, 2002; Thomson and Bannister, 2003; Binzegger et al., 2004). A more realistic graph structure might eliminate the “fan-out” tendency that underlies the instability of the fixed point. Such adaptations may also entail the use of a more complex neuron model than assumed here. Kumar et al. (2008) demonstrated that the regime enabling stable propagation activity in a feed-forward sub-network embedded in a locally connected random network is much greater for a conductance-based neuron model than previously found for a current-based model (Mehring et al., 2003). Moreover, a more sophisticated neuron model, particularly with respect to the dendritic integration of inputs, might exhibit greater symmetry breaking properties.

Another candidate is the formulation of the STDP rule. Although we have shown that our results are robust to the choice of update rule, we have only considered extremely simplified models of STDP exclusively at excitatory–excitatory synapses; the role of plasticity at inhibitory synapses is neglected. We addressed models in which the STDP window is fitted by two exponential functions, which is easy to analyze but may be too simplistic. Further, we have so far only considered rules based on pairs of spikes. Pfister and Gerstner (2006) introduced a rule based on triplets of spikes which accounts for experimental findings on the dependence of the weight change on the frequency of the pairing protocols (Sjöström et al., 2001). The triplet model has been shown to map to the BCM rule (Bienenstock et al., 1982) and to give near-optimal information transmission when combined with spike-frequency adaptation in the post-synaptic neuron (Hennequin et al., 2010). The triplet interaction could well also prove to be important for the development of structure. Moreover, recent experimental findings reveal that neuromodulation can influence STDP. Seol et al. (2007) showed that the polarity of the STDP window depends on the activation of receptors sensitive to cholinergic and adrenergic activity. A recent study by Kozloski and Cecchi (2010) demonstrated that polysynaptic loops eliminated by standard STDP can be restored by such a reversed rule. This finding suggests that the interaction of causal and anti-causal windows may indeed result in the kind of symmetry breaking that would enable groups of neurons to develop a qualitatively different response to a given synchronous stimulus. Neuromodulation can also be essential either for inducing STDP or for altering the threshold for induction, see Pawlak et al. (2010) in this special issue for a review. Finally, our theoretical approach does not take into account interactions of STDP with other forms of plasticity except post-synaptic homeostasis. Although no development of feed-forward structure occurs in our homeostatically regulated network models, the presence of homeostasis certainly has a strong stabilizing effect on the network dynamics. Therefore, homeostasis may be a crucial element that enables network dynamics to be maintained during stimulus-driven development of structure. Short-term plasticity is another feature overlooked in the current study. The interaction of STDP with depression and facilitation could result in different network dynamics that would be more favorable to symmetry breaking, as could the incorporation of tagging processes to stabilize the developed structures (Frey and Morris, 1997; Reymann and Frey, 2007; Clopath et al., 2008). We therefore conclude that the development and investigation of more sophisticated models of STDP that interact with other forms of plasticity on different time scales would be a useful future extension of this research.

A third area in which the model presented here could be adapted is the stimulus protocol. Here, we investigated whether structure can develop from the seed of a single correlated group. Several forms of synaptic pathway organization have been identified and mathematically analyzed for recurrent networks under the assumption that they receive input from more than one pool of correlated inputs (e.g., Gilson et al., 2009b,c; see Gilson et al., 2010 for a thorough review). However, the majority of the findings assume STDP that is independent of the synaptic weight or only very weakly dependent. Moreover, only small networks are considered in which each neuron receives external correlated input. It will therefore be interesting to discover how these results scale up to the large, sparsely connected networks with weight-dependent plasticity considered in this article.

These candidate adaptations are not intended to be exhaustive. There are many ways in which the theoretical approach and network model presented here could be improved and supplemented. However, we hope that the suggestions above can serve to stimulate research to reveal the necessary and sufficient conditions for the development of structure in the mammalian brain.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

Partially funded by DIP F1.2, BMBF Grant 01GQ0420 to BCCN Freiburg, EU Grant 15879 (FACETS), EU Grant 269921 (BrainScaleS), Helmholtz Alliance on Systems Biology (Germany), Next-Generation Supercomputer Project of MEXT (Japan), Neurex, and the Junior Professor Program of Baden-Württemberg. High-performance computing facilities were made available through the Norwegian Metacenter for Computational Science (Notur) and JUGENE Grant JINB33 at the Research Center Jülich. We gratefully acknowledge T. Potjans for his help in determining appropriate connectivity levels for our network model and S. Jahnke for his comments on an earlier version of the manuscript.

## References

Abbott, L. F., and Nelson, S. B. (2000). Synaptic plasticity: taming the beast. *Nat. Neurosci. *3(Suppl.), 1178–1183.

Abeles, M. (1991). *Corticonics: Neural Circuits of the Cerebral Cortex*, 1st Edn. Cambridge: Cambridge University Press.

Abeles, M., Bergman, H., Margalit, E., and Vaadia, E. (1993). Spatiotemporal firing patterns in the frontal cortex of behaving monkeys. *J. Neurophysiol. *70, 1629–1638.

Bi, G.-Q., and Poo, M.-M. (1998). Synaptic modifications in cultured hippocampal neurons: dependence on spike timing, synaptic strength, and postsynaptic cell type. *J. Neurosci. *18, 10464–10472.

Bienenstock, E. L., Cooper, L. N., and Munro, P. W. (1982). Theory for the development of neuron selectivity: orientation specificity and binocular interaction in visual cortex. *J. Neurosci. *2, 32–48.

Binzegger, T., Douglas, R. J., and Martin, K. A. C. (2004). A quantitative map of the circuit of cat primary visual cortex. *J. Neurosci. *39, 8441–8453.

Brunel, N. (2000). Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. *J. Comput. Neurosci. *8, 183–208.

Clopath, C., Ziegler, L., Vasilaki, E., Büsing, L., and Gerstner, W. (2008). Tag-trigger-consolidation: a model of early and late long-term-potentiation and depression. *PLoS Comput. Biol. *4, e1000248. doi: 10.1371/journal.pcbi.1000248

Debanne, D., Gähwiler, B. H., and Thompson, S. M. (1996). Cooperative interactions in the induction of long-term potentiation and depression of synaptic excitation between hippocampal CA3-CA1 cell pairs in vitro. *Proc. Natl. Acad. Sci. U.S.A. *93, 11225–11230.

Debanne, D., Gähwiler, B. H., and Thompson, S. M. (1999). Heterogeneity of synaptic plasticity at unitary CA1-CA3 and CA3-CA3 connections in rat hippocampal slice cultures. *J. Neurosci. *19, 10664–10671.

Diesmann, M., Gewaltig, M.-O., and Aertsen, A. (1999). Stable propagation of synchronous spiking in cortical neural networks. *Nature *402, 529–533.

Doursat, R., and Bienenstock, E. (2006). “Neocortical self-structuration as a basis for learning,” in *5th International Conference on Development and Learning (ICDL 2006)*, Indiana University, Bloomington, IN.

Eckhorn, R., Bauer, R., Jordan, W., Brosch, M., Kruse, W., Munk, M., and Reitböck, H. J. (1988). Coherent oscillations: a mechanism of feature linking in the visual cortex? *Biol.* *Cybern. *60, 121–130.

Fiete, I. R., Senn, W., Wang, C. Z. H., and Hahnloser, R. H. R. (2010). Spike-time-dependent plasticity and heterosynaptic competition organize networks to produce long scale-free sequences of neural activity. *Neuron *65, 563–576.

Gerstner, W., Ritz, R., and van Hemmen, J. (1993). Why spikes? Hebbian learning and retrieval of time-resolved excitation patterns. *Biol. Cybern. *69, 503–515.

Gewaltig, M.-O., Diesmann, M., and Aertsen, A. (2001). Propagation of cortical synfire activity: survival probability in single trials and stability in the mean. *Neural Netw. *657–673.

Gilson, M., Burkitt, A. N., Grayden, D. B., Thomas, D. A., and van Hemmen, J. L. (2009a). Emergence of network structure due to spike-timing-dependent plasticity in recurrent neuronal networks. I. Input selectivity – strengthening correlated input pathways.* Biol. Cybern. *101, 81–102.

Gilson, M., Burkitt, A. N., Grayden, D. B., Thomas, D. A., and van Hemmen, J. L. (2009b). Emergence of network structure due to spike-timing-dependent plasticity in recurrent neuronal networks. II. Input selectivity – symmetry breaking.* Biol. Cybern. *101, 103–114.

Gilson, M., Burkitt, A. N., Grayden, D. B., Thomas, D. A., and van Hemmen, J. L. (2009c). Emergence of network structure due to spike-timing-dependent plasticity in recurrent neuronal networks IV. *Structuring synaptic pathways among recurrent connections. Biol. Cybern. *101, 427–444.

Gilson, M., Burkitt, A. N., and van Hemmen, J. L. (2010). STDP in recurrent neuronal networks. *Front. Comput. Neurosci. *4:23. doi: 10.3389/fncom.2010.00023

Gray, C. M., and Singer, W. (1989). Stimulus-specific neuronal oscillations in orientation columns of cat visual cortex. *Proc. Natl. Acad. Sci. U.S.A. *86, 1698–1702.

Gustafsson, B., Wigstrom, H., Abraham, W. C., and Huang, Y.-Y. (1987). Long-term potentiation in the hippocampus using depolarizing current pulses as the conditioning stimulus to single volley synaptic potentials. *J. Neurosci. *7, 774–780.

Gütig, R., Aharonov, R., Rotter, S., and Sompolinsky, H. (2003). Learning input correlations through nonlinear temporally asymmetric Hebbian plasticity. *J. Neurosci. *23, 3697–3714.

Hahnloser, R. H., Kozhevnikov, A. A., and Fee, M. S. (2002). An ultra-sparse code underlies the generation of neural sequences in a songbird. *Nature *419, 65–70.

Hayon, G., Abeles, M., and Lehmann, D. (2004). Modeling compositionality by dynamic binding of synfire chains. *J. Comput. Neurosci. *17, 179–201.

Hebb, D. O. (1949). *The Organization of Behavior: A Neuropsychological Theory*. New York: John Wiley and Sons.

Hennequin, G., Gerstner, W., and Pfister, J.-P. (2010). STDP in adaptive neurons gives close-to-optimal information transmission. *Front. Comput. Neurosci. *4:143. doi: 10.3389/fncom.2010.00143

Hertz, J., and Prügel-Bennet, A. (1996). Learning synfire-chains by self-organization. *Network *7, 357–363.

Hosaka, R., Araki, O., and Ikeguchi, T. (2008). STDP provides the substrate for igniting synfire chains by spatiotemporal input patterns. *Neural Comput. *20, 415–435.

Iglesias, J., Eriksson, J., Grize, F., Tomassini, M., and Villa, A. (2005). Dynamics of pruning in simulated large-scale spiking neural networks. *Biosystems *79, 11–20.

Ikegaya, Y., Aaron, G., Cossart, R., Aronov, D., Lampl, I., Ferster, D., and Yuste, R. (2004). Synfire chains and cortical songs: temporal modules of cortical activity. *Science *5670, 559–564.

Izhikevich, E. M., Gally, J. A., and Edelman, G. M. (2004). Spike-timing dynamics of neuronal groups. *Cereb. Cortex *14, 933–944.

Jun, J. K., and Jin, D. Z. (2007). Development of neural circuitry for precise temporal sequences through spontaneous activity, axon remodeling, and synaptic plasticity. *PLoS ONE *2, e723. doi: 10.1371/journal.pone.0000723

Kempter, R., Gerstner, W., and van Hemmen, J. L. (1999). Hebbian learning and spiking neurons. *Phys. Rev. E *59, 4498–4514.

Kempter, R., Gerstner, W., and van Hemmen, J. L. (2001). Intrinsic stabilization of output rates by spike-based Hebbian learning. *Neural Comput. *12, 2709–2742.

Kistler, W. M., and van Hemmen, J. L. (2000). Modeling synaptic plasticity in conjunction with the timing of pre- and postsynaptic action potentials. *Neural Comput. *12, 385–405.

Kozhevnikov, A., and Fee, M. S. (2007). Singing-related activity of identified HVC neurons in the zebra finch. *J. Neurophysiol. *97, 4271–4283.

Kozloski, J., and Cecchi, G. A. (2010). A theory of loop formation and elimination by spike timing-dependent plasticity. *Front. Neural Circuits *4:7. doi: 10.3389/fncir.2010.00007

Kuhn, A., Aertsen, A., and Rotter, S. (2003). Higher-order statistics of input ensembles and the response of simple model neurons. *Neural Comput. *1, 67–101.

Kumar, A., Rotter, S., and Aertsen, A. (2008). Conditions for propagating synchronous spiking and asynchronous firing rates in a cortical network model. *J. Neurosci. *28, 5268–5280.

Levy, N., Horn, D., Meilijson, I., and Ruppin, E. (2001). Distributed synchrony in a cell assembly of spiking neurons. *Neural Netw. *14, 815–824.

Lewis, D., Melchitzky, D., and Burgos, G.-G. (2002). Specificity in the functional architecture of primate prefontal cortex. *J. Neurocytol. *31, 265–276.

Liu, J. K., and Buonomano, D. V. (2009). Embedding multiple trajectories in simulated recurrent neural networks in a self-organizing manner. *J. Neurosci. *29, 13172–13181.

Lubenov, E. V., and Siapas, A. G. (2008). Decoupling through synchrony in neuronal circuits with propagation delays. *Neuron *58, 118–131.

Markram, H., Lübke, J., Frotscher, M., and Sakmann, B. (1997). Regulation of synaptic efficacy by coincidence of postsynaptic APs and EPSPs. *Science *275, 213–215.

Markram, H., and Sakmann, B. (1995). Action potentials propagating back into dendrites trigger changes in efficacy of single-axon synapses between layer V pyramidal neurons. *Soc. Neurosci. Abstr. *21, 2007.

Masuda, N., and Kori, H. (2007). Formation of feedforward networks and frequency synchrony by spike-timing-dependent plasticity. *J. Comput. Neurosci. *22, 327–345.

Mehring, C., Hehl, U., Kubo, M., Diesmann, M., and Aertsen, A. (2003). Activity dynamics and propagation of synchronous spiking in locally connected random networks. *Biol. Cybern. *88, 395–408.

Migliore, M., Cannia, C., Lytton, W. W., Markram, H., and Hines, M. (2006). Parallel network simulations with neuron. *J. Comput. Neurosci. *21, 119–223.

Montgomery, J., Pavlidis, P., and Madison, D. (2001). Pair recordings reveal all-silent synaptic connections and the postsynaptic expression of long-term potentiation. *Neuron *29, 691–701.

Morrison, A., Aertsen, A., and Diesmann, M. (2007). Spike-timing dependent plasticity in balanced random networks. *Neural Comput. *19, 1437–1467.

Morrison, A., Diesmann, M., and Gerstner, W. (2008). Phenomenological models of synaptic plasticity based on spike-timing. *Biol. Cybern. *98, 459–478.

Morrison, A., Mehring, C., Geisel, T., Aertsen, A., and Diesmann, M. (2005). Advancing the boundaries of high connectivity network simulation with distributed computing. *Neural Comput. *17, 1776–1801.

Nordlie, E., Plesser, H. E., and Gewaltig, M.-O. (2009). Towards reproducible descriptions of neuronal network models. *PLoS Comput. Biol. *5, e1000456. doi: 10.1371/journal.pcbi.1000456

Pawlak, V., Wickens, J. R., Kirkwood, A., and Kerr, J. N. D. (2010). Timing is not everything: neuromodulation opens the STDP gate. *Front. Syn. Neurosci. *2:146. doi: 10.3389/fnsyn.2010.00146

Pecevski, D., Natschläger, T., and Schuch, K. (2009). PCSIM: a parallel simulation environment for neural circuits fully integrated with python. *Front. Neuroinform. *3:11. doi: 10.3389/neuro.11.011.2009

Pfister, J.-P., and Gerstner, W. (2006). Triplets of spikes in a model of spike timing-dependent plasticity. *J. Neurosci. *26, 9673–9682.

Prut, Y., Vaadia, E., Bergman, H., Haalman, I., Hamutal, S., and Abeles, M. (1998). Spatiotemporal structure of cortical activity: Properties and behavioral relevance. *J. Neurophysiol. *79, 2857–2874.

Pulvermüller, F., and Shtyrov, Y. (2009). Spatiotemporal signatures of large-scale synfire chains for speech processing as revealed by MEG. *Cereb. Cortex *19, 79–88.

Reymann, K. G., and Frey, J. U. (2007). The late maintenance of hippocampal LTP: Requirements, phases, “synaptic tagging,” “late-associativity” and implications. *Neuropharmacology *52, 24–40.

Rubin, J., Lee, D., and Sompolinsky, H. (2001). Equilibrium properties of temporally asymmetric Hebbian plasticity. *Phys. Rev. Lett. *86, 364–367.

Schüz, A., and Braitenberg, V. (2002). “The human cortical white matter: quantitative aspects of cortico-cortical long-range connectivity,” in *Cortical Areas: Unity and Diversity*, Chapter 16, eds A. Schüz and R. Miller (London: Taylor and Francis), 377–385.

Seol, G., Ziburkus, J., Huang, S., Song, L., Kim, I., Takamiya, K., Huganir, R., Lee, H.-K., and Kirkwood, A. (2007). Neuromodulators control the polarity of spike-timing-dependent plasticity. *Neuron *55, 919–929.

Sjöström, P., Turrigiano, G., and Nelson, S. (2001). Rate, timing, and cooperativity jointly determine cortical synaptic plasticity. *Neuron *32, 1149–1164.

Song, S., Miller, K. D., and Abbott, L. F. (2000). Competitive Hebbian learning through spike-timing-dependent synaptic plasticity. *Nat. Neurosci. *3, 919–926.

Song, S., Sjöström, P. J., Reigl, M., Nelson, S., and Chklovskii, D. B. (2005). Highly nonrandom features of synaptic connectivity in local cortical circuits. *PLoS Biol. *3, e350. doi: 10.1371/journal.pbio.0030350

Standage, D., Jalil, S., and Trappenberg, T. (2007). Computational consequences of experimentally derived spike-time and weight dependent plasticity rules. *Biol. Cybern. *96, 615–623.

Thomson, A. M., and Bannister, P. (2003). Interlaminar connections in the neocortex. *Cereb. Cortex *13, 5–14.

Turrigiano, G. G., Leslie, K. R., Desai, N. S., Rutherford, L. C., and Nelson, S. B. (1998). Activity-dependent scaling of quantal amplitude in neocortical neurons. *Nature *391, 892–896.

van Rossum, M. C. W., Bi, G.-Q., and Turrigiano, G. G. (2000). Stable Hebbian learning from spike timing-dependent plasticity. *J. Neurosci. *20, 8812–8821.

Wang, H.-X., Gerkin, R. C., Nauen, D. W., and Bi, G.-Q. (2005). Coactivation and timing-dependent integration of synaptic potentiation and depression. *Nat. Neurosci. *8, 187–193.

Keywords: spike-timing dependent plasticity, feed-forward structure, random recurrent network

Citation: Kunkel S, Diesmann M and Morrison A (2011) Limits to the development of feed-forward structures in large recurrent neuronal networks. *Front. Comput. Neurosci.* **4**:160. doi: 10.3389/fncom.2010.00160

Received: 15 March 2010;
Accepted: 25 December 2010;

Published online: 14 February 2011.

Edited by:

Per Jesper Sjöström, University College London, UKReviewed by:

Brent Doiron, University of Pittsburgh, USASean L. Hill, Ecole Polytechnique Fédérale de Lausanne, Switzerland

Tim P. Vogels, Ecole Polytechnique Fédérale de Lausanne, Switzerland

Copyright: © 2011 Kunkel, Diesmann and Morrison. This is an open-access article subject to an exclusive license agreement between the authors and Frontiers Media SA, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are credited.

*Correspondence: Abigail Morrison, Bernstein Center Freiburg, Hansastrasse 9a, 79104 Freiburg im Breisgau, Germany. e-mail: morrison@bcf.uni-freiburg.de