{'fraction_non_alphanumeric': 0.03368987015958359, 'fraction_numerical': 0.01630155007721787, 'mean_word_length': 4.402812113720643, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Self-organization of feedforward structure and entrainment in excitatory neural networks with spike-timing-dependent plasticity

Spike-timing dependent plasticity (STDP) is an organizing principle of biological neural networks. While synchronous firing of neurons is considered to be an important functional block in the brain, how STDP shapes neural networks possibly toward synchrony is not entirely clear. We examine relations between STDP and synchronous firing in spontaneously firing neural populations. Using coupled heterogeneous phase oscillators placed on initial networks, we show numerically that STDP prunes some synapses and promotes formation of a feedforward network. Eventually a pacemaker, which is the neuron with the fastest inherent frequency in our numerical simulations, emerges at the root of the feedforward network. In each oscillatory cycle, a packet of neural activity is propagated from the pacemaker to downstream neurons along layers of the feedforward network. This event occurs above a clear-cut threshold value of the initial synaptic weight. Below the threshold, neurons are self-organized into separate clusters each of which is a feedforward network.


I. INTRODUCTION
Synchronous firing of neurons has been widely observed and is considered to be a neural code that adds to firing rates. For example, experimental evidence suggests the relevance of synchronous firing in stimulus encoding [1], feature binding [2,3], and selective attention [3,4]. Collective dynamical states of neurons including synchrony may appear as a result of self-organization based on synaptic plasticity. Modification of synaptic weights (i.e., weights of edges in the network terminology) often occurs in a manner sensitive to relative spike timing of presynaptic and postsynaptic neurons, which is called spike-timing-dependent plasticity (STDP). In the commonly found asymmetric STDP, which we consider in this work, long-term potentiation (LTP) occurs when presynaptic firing precedes postsynaptic firing by tens of milliseconds or less, and long-term depression (LTD) occurs in the opposite case [5]. The amount of plasticity is larger when the difference in the presynaptic spike time and the postsynaptic spike time is smaller [5].
The asymmetric STDP reinforces causal pairs of presynaptic and postsynaptic spikes and eliminate other pairs. Based on this property of STDP, how STDP may lead to various forms of synchronous firing has been studied in both experiments and theory. Synchronous firing in the sense of simultaneity of spike timing can be established in recurrent neural networks when the strength of LTP and that of LTD are nearly balanced [6]. Large-scale numerical simulations suggest that reproducible spatiotemporal patterns of spike trains self-organize in heterogeneous recurrent neural networks [7,8]. Self-organization of clusters of synchronously firing neurons that excite each other in a cyclic manner has also been reported [9,10].
We previously showed that STDP leads to formation of feedforward networks and entrainment when there is a pacemaker in the initial network [11]. We considered random networks of coupled oscillators whose synaptic weights change slowly via STDP. We assumed that the oscillators have a common inherent frequency except a single pacemaker whose inherent frequency is larger. By definition, the rhythm of the pacemaker is not affected by those of other oscillators. The network generated via STDP is a feedforward network whose root is the pacemaker. In a final network, a spike packet travels from the pacemaker to the other neurons in a laminar manner. The neurons directly postsynaptic to the pacemaker fire more or less synchronously just after the pacemaker does. These neurons form the first layer. These neurons induce synchronous firing of the neurons directly postsynaptic to them, which define the second layer. In this fashion, a spike packet starting from the pacemaker reaches the most downstream neurons within relatively short time, which resembles the phenomenology of the synfire chain [12]. Compared to the case of frozen synaptic weights, a pacemaker entrains the rest of the network more easily with STDP in the meaning that entrainment occurs with smaller initial synaptic weights.
The previous work does not explain how pacemakers emerge. No matter whether the pacemakers are intrinsic oscillators or network oscillators, they pace rhythms of other elements without being crucially affected by other rhythms. Although some pacemakers may be "robust" oscillators whose rhythms are insensitive to general input, a more natural explanation may be that pacemakers emerge through synaptic plasticity in a neural network in which pacemakers are initially absent. In this case, emergent pacemakers do not have to be robust oscillators; their rhythms can change in response to external input. The emergent network topology makes such neurons pacemakers by eliminating incoming synapses. A neuron would fire with its own rhythm if it is not downstream to any neuron. This scenario is actually the case for two-neuron networks [11,13]. Here we are concerned to networks of more than two neurons. An associated question is which oscillator may become a pacemaker.
In this work, we numerically investigate recurrent networks of coupled phase oscillators subject to STDP. We show that, when the initial synaptic weights are strong enough, STDP indeed yields feedforward networks so that downstream neurons are entrained by an emergent pacemaker. To our numerical evidence, the emergent pacemaker is always the neuron with the largest intrinsic frequency. Below the threshold for entrainment, STDP leads to the segregation of the initial neural network into subnetworks of feedforward networks.

A. Coupled phase oscillators
We model dynamics of neural networks by N coupled phase oscillators whose synaptic weights are plastic. Although a majority of real neurons fire in the excitable (i.e., fluctuationdriven) regime, for tractability we use phase oscillators, which fire in an oscillatory manner.
Generally speaking, phase transitions are more easily and clearly determined in the oscillatory regime than in the excitable regime. This is a reason why collective neural dynamics [14,15] including ones associated with STDP [6,11,16] have been analyzed in the oscillatory regime actually to give insights into dynamics of neural networks possibly operating in the excitable regime. In the following, we report numerical results for N = 3 and N = 100.
We identify φ i = 0 and φ i = 2π. When φ i crosses 0 in the positive direction, neuron i is defined to fire. We denote by t j and t i the spike time of presynaptic and postsynaptic neurons. If φ i crosses 0 in the positive direction as time advances from t to t + ∆t, we set As the initial condition, we set φ i = 0 (1 ≤ i ≤ N) for N = 3. We adopt this artificial initial condition to draw phase diagrams to systematically understand possible routes to synchrony via STDP. For N = 100, φ i (0) is picked randomly and independently for each i from the uniform density on [0, 2π). Neuron i is endowed with inherent frequency ω i so that it fires regularly at rate ω i /2π when isolated.
Connectivity between neurons is unidirectional and weighted, consistent with the properties of chemical synapses. The set of edges in a network is denoted by E. In other words, (j, i) ∈ E if neuron j is presynaptic to neuron i. Dynamics of the coupled phase oscillators are given by where k is the average in-degree of neuron i, g ji is a synaptic weight, and ξ i represents the standard Gaussian white noise independent for different i. As a result of the phase reduction theory [17], the coupling term in the oscillatory regime is generally given by a 2π-periodic function of the phase difference φ j − φ i under the assumption of weak coupling. This is also the case for pulse coupling, for which averaging an original pulse coupling term over one oscillatory cycle results in a coupling term as a function of φ j − φ i [14]. Modeling realistic synaptic coupling needs a coupling term that contains higher harmonics [14]. However, our objective in the present paper is not to precisely describe the neural dynamics but to clarify general consequences of STDP under the oscillatory condition. We thus employ the simplest coupling term (i.e sinusoidal coupling).
For N = 3, we set the amplitude of the noise σ = 0.0071 so that the phase transitions are sharp enough and artificial resonance that is prone to occur when inherent frequencies satisfy M i ω i = M j ω j for small integers M i and M j (1 ≤ i < j ≤ N) is avoided. Accordingly, an independent normal variable with mean 0 and standard deviation σ √ ∆t = 0.00071 is added to each neuron every time step ∆t = 0.01; we use the Euler-Maruyama integration scheme with unit time ∆t. To determine the phase transitions for N = 100, we do not apply dynamical noise because, up to our numerical efforts, the numerical results do not significantly suffer from artificial resonance. In some other simulations with N = 100, we add different amplitudes of dynamical noise to examine the robustness of the results.

B. STDP
With STDP, g ji is repeatedly updated depending on spike timing of neuron j and i.
Specifically, LTP occurs when a postsynaptic neuron fires slightly after a presynaptic neuron does, and LTD occurs in the opposite case [5]. We assume that synaptic plasticity operates much more slowly than firing dynamics. We denote by A + and A − the maximum amount of LTP and that of LTD incurred by a single STDP event. Most of previous theoretical work supposes that A − is somewhat, but not too much, larger than A + , to avoid explosion in firing rates and to keep neurons firing [6,7,8,9,10,11]. Therefore we set A + /A − = 0.9. How a single spike pair specifically modifies the synaptic weight is under investigation [8,18], and triplets or higher-order combination of presynaptic and postsynaptic spikes rather than a single presynaptic and postsynaptic spike pair may induce STDP [19]. However, we consider the simplest situation in which STDP modifies synaptic weights in an additive manner and the amount of STDP is determined by the relative timing of a presynaptic and postsynaptic spike pair. A single synaptic modification ∆g ji triggered by a spike pair is represented by where τ is the characteristic time scale of the learning window, which is known in experiments to be 10-20 ms [5]. Given that inherent frequencies of many pyramidal neurons roughly range between 5 and 20 Hz, τ is several times smaller than a characteristic average interspike interval. Therefore, following [11], we set τ = 1/6 × 2π/ω, where ω is a typical value of spike frequency that is used to determine ω i . Following our previous work [11], we set ω = 8.1.
Because learning is slow compared to neural dynamics, A − must be by far smaller than a typical value of g. To satisfy this condition, we set A − = 0.001 for N = 3. When N = 100, average in-degree k is set equal to 10. This implies that a neuron receives about five to ten times more synapses than when N = 3. To normalize this factor, we set it A − = 0.0001 for N = 100.
We assume that g ji is confined in [0, g max ]; all the synapses are assumed to be excitatory, because the asymmetric STDP explained in Sec. I has been found mostly in excitatory synapses. Because dynamical noise is assumed not to be large, all the synaptic weights usually develop until g ji almost reaches either g max or 0, until when we run each simulation run. Note that, even if g ji = 0 is reached, (j, i) still belongs to E. The upper limit g max is determined so that a notion of synchronization that we define in Sec. II C does not occur when g ji = g max , ∀(j, i) ∈ E. Accordingly, we set g max = 7.5 and g max = 15 for N = 3 and N = 100, respectively.

C. Measurement of synchrony
To obtain the threshold for synchrony in Secs. III A and III B 1, we start numerical simulations with the initial condition g ji = g 0 , ∀(j, i) ∈ E. There are various notions of synchrony.
We focus on the possibility of frequency synchrony in which neurons fire at the same rate.
In the oscillatory regime, frequency synchrony is commonly achieved in two main ways. One is when neurons are connected by sufficiently strong mutual coupling. Then they oscillate at the same rate and with proximate phases. The other is when some neurons entrain others.
When upstream neurons, which serve as pacemakers, entrain downstream neurons so that they are synchronized in frequency, synchronous firing in the sense of spike timing may be missing due to synaptic delay. However, neurons located at the same level in the hierarchy relative to the pacemakers tend to have close spike timing [11,20]. We explore possible emergence of such dynamics when pacemakers are initially absent in networks.
We quantify the degree of frequency synchrony by order parameter r defined by whereω i = dφ i /dt is the actual instantaneous frequency of neuron i when coupled to other neurons. If all the neurons fire exactly at the same rate, r would become negative infinity.
In the actual frequency synchrony, r takes a large negative value mainly because of time discretization. We manually set r c = −4 for N = 3 and r c = −9 for N = 100, so that r ≤ r c corresponds to the full frequency synchrony. The value of r c for N = 100 is smaller than for N = 3 for two reasons. First, in the numerical simulations determining the degree of frequency synchrony, dynamical noise is present for N = 3 and absent for N = 100. Second, we are concerned to the frequency synchrony of all the neurons so that i (ω i −1/N i ′ωi ′ ) 2 is small regardless of N; we have to normalize the prefactor 1/N in Eq. (3).

A. Networks of three neurons
Our goal is to understand dynamics of large neural networks. As a starting point, we examine network evolution and possibility of frequency synchrony using small networks, which will help us understand dynamics of large networks. Two-neuron networks were previously analyzed [11]. We need at least three neurons to understand competition between different synapses, pruning of synapses, and effects of heterogeneity. Accordingly, we examine dynamics of different three-neuron networks under STDP.

Complete graph
Consider the complete graph [ Fig. 1(a)], in which every pair of neurons is bidirectionally connected. The complete graph does not survive STDP because LTP of a synapse implies LTD of the synapse in the reversed direction and the amount of LTD is assumed to be larger than that of LTP for the same time lag. We examine which synapses survive and whether frequency synchrony emerges through STDP. If a predetermined pacemaker exists in a network, the activity of the other neurons will be entrained into the rhythm of the pacemaker with sufficiently large initial synaptic weights, which was previously shown for N = 2 and N = 100 [11]. Here we consider N = 3 and compare numerical results when a pacemaker is initially present and absent in the complete graph. Note that the effective initial network when the pacemaker neuron 1 is initially present is the one shown in Fig. 1(b), because the synapses toward the pacemaker are defined to be entirely ineffective.
First, we examine the relation between heterogeneity in inherent frequencies, initial synaptic weights, and synchrony. We expect that small heterogeneity and large initial synaptic weights favor synchrony. To focus on phase transitions, we reduce the number of parameters by setting all the initial synaptic weights equal to g 0 and restrain inherent frequencies ω 1 , ω 2 , and ω 3 (ω 1 ≥ ω 2 ≥ ω 3 ) by imposing, ω 1 − ω 2 = ω 2 − ω 3 ≡ ∆ω, where ω 2 = 8.1. Numerically obtained phase diagrams are shown in Figs. 2(a) and 2(b) for the cases in which a pacemaker is initially present and absent, respectively. The results are qualitatively the same for the two situations. The neurons get disconnected and fire independently as a result of STDP for sufficiently small g 0 or sufficiently large ∆ω (blue regions labeled D). A feedforward network whose root is the fastest oscillator emerges for sufficiently large g 0 or sufficiently small ∆ω (yellow regions labeled A). Then all the neurons rotate at frequency ω 1 . In the intermediate regime (green regions labeled C), final synaptic weights satisfy g 23 ≈ g max and g 12 , g 13 , g 21 , g 31 , and g 32 ≈ 0. In this case, neuron 2 entrains neuron 3 so that they oscillate at frequency ω 2 , whereas neuron 1 gets disconnected and oscillates at frequency ω 1 . We rarely observed the case in which neuron 1 entrains 2 (or 3) and neuron 3 (or 2) gets isolated. Although ω 1 − ω 2 = ω 2 − ω 3 , neuron 1 is more likely to segregate from the network than neuron 3 is. Quantitatively speaking, Figs. 2(a) and 2(b) indicate that the entrainment of the entire network by the fastest neuron (i.e., neuron 1) is to some extent easier to realize when the pacemaker is initially absent than present (yellow regions labeled A). In Figs. 2(a) and 2(b), the phase diagrams are disturbed along vertical lines at ∆ω ≈ 2.7. This artifact comes from the fact that ω 1 , ω 2 , and ω 3 approximately satisfy the resonance condition (i.e., M 1 ω 1 = M 2 ω 2 = M 3 ω 3 with small integers M 1 , M 2 , and M 3 ). In some of the following figures, similar disturbance appears along special lines. We can wash away these artifacts by increasing the amount of dynamical noise. However, we prefer not doing so to prevent the boundaries between different phases from being blurred too much.
Next, to examine what happens when ω 1 , ω 2 , and ω 3 change independently, we set g 0 = 0.15, ω 2 = 8.1, and vary ∆ω 1 ≡ ω 1 − ω 2 and ∆ω 2 ≡ ω 2 − ω 3 . Numerical results with and To examine the cause of the asymmetry in Figs. 2(c) and 2(d) along the two lines ∆ω 1 = 0 and ∆ω 2 = 0, we analyze a two-neuron network with asymmetric initial synaptic weights shown in Fig. 3(a). The two neurons h and l have inherent frequency ω h and ω l (≤ ω h ). The weights of the synapse from neuron h to neuron l and that from neuron l to neuron h are denoted by g f and g b , respectively. When ∆ω 1 = 0 and ∆ω 2 ≥ 0 in the three-neuron network, neurons 1 and 2 are synchronized almost from the beginning, in both frequency and phase, because ω 1 = ω 2 . This is true if a trivial condition g 12 +g 21 > 0 is satisfied. Then the network is reduced to the two-neuron network by identifying ω h = ω 1 = ω 2 , ω l = ω 3 , g f = g 13 + g 23 , and g b = (g 21 + g 31 )/2. When, ∆ω 1 ≥ 0 and ∆ω 2 = 0 in the three-neuron network, neurons 2 and 3 are synchronized in frequency and phase as far as g 23 + g 32 > 0. Then the network is reduced to the two-neuron network with ω h = ω 1 , ω l = ω 2 = ω 3 , g f = (g 12 + g 13 )/2, and g b = g 21 + g 31 . For these two situations, we calculate the threshold for frequency synchrony in the two-neuron network using the semi-analytical method developed in [11]. Because all the synaptic weights are initially equal to g 0 in Fig. 2, the initial condition for the two- for ∆ω 1 ≡ ∆ω ≥ 0, ∆ω 2 = 0. The phase-transition curves for the frequency synchrony are shown in Fig. 3(b), indicating that the threshold is larger along the ∆ω 2 = 0 line than along the ∆ω 1 = 0 line. This is consistent with the three-neuron results shown in Figs. 2(c) and 2(d).

Feedforward loop
Other three-neuron networks, particularly feedforward ones, are presumably embedded in larger neural networks in the course of network evolution. First, we consider the network shown Fig. 4(a) as the initial network. The phase diagram in the ∆ω 1 -∆ω 2 parameter space with g 0 = 0.15 is shown in Fig. 4(b), which looks similar to Figs. 2(c) and 2(d). As in the case of the complete graph, the situation in which neuron 1 entrains neuron 2 with neuron 3 isolated is less likely to arise than that in which neuron 2 entrains neuron 3 with neuron 1 isolated.

Fan-in network
Next, we examine dynamics starting from the fan-in network shown in Fig. 5(a). In this network, neuron 3 is postsynaptic to two pacemaker neurons 1 and 2. We are concerned to which neuron entrains neuron 3.
First, we examine the case in which two synapses are initially equally strong and the inherent frequencies of the two upstream neurons are different. Accordingly we set g 13 = when ω 1 ≈ ω 2 (yellow, A). This phase extends to the disconnection phase (blue, D) on the diagonal because, on this line, the firing of neuron 1 elicits LTP of both synapses so does firing of neuron 2. However, this situation is not generic in that ω 1 and ω 2 must be very close for this to happen. When ω 1 and ω 2 are not close to each other and not too far from ω 3 , which upstream neuron entrains neuron 3 is not obvious. Figure 5(b) tells that a necessary condition for an upstream neuron to entrain neuron 3 is that the difference between its inherent frequency and ω 3 is less than ≈ 1.0. This condition roughly corresponds to the requirement for the entrainment in the two-neuron feedforward network with g 0 = 0.2. This explains the two rectangular regions ∆ω 1 > 1.0, ∆ω 2 < 1.0, and ∆ω 1 < 1.0, ∆ω 2 > 1.0 of Fig. 5(b). In the remaining region (i.e. ∆ω 1 < 1.0 and ∆ω 2 < 1.0), the upstream neuron whose inherent frequency is closer to ω 3 , equivalently, the slower upstream neuron, largely wins the competition (regions marked by ). The faster upstream neuron entrains neuron 3 when the inherent frequency of the slower upstream neuron is very close to ω 3 (regions marked by ). The total size of the latter regions is much smaller than that of the former regions.

B. Networks of many neurons
In this section, we use networks of heterogeneous N = 100 neurons to examine what network structure and dynamics self-organize via STDP when we start from random neural networks. The inherent frequencies of the neurons are independently picked from the truncated Gaussian distribution with mean 8.1, standard deviation 0.5, and support ω i ∈ [7.6, 8.6].
We assume that every neuron has k = 10 randomly selected presynaptic neurons on average so that an arbitrary pair of neurons is connected by a directed edge with probability k /(N − 1) ≈ 0.1. Except in Sec. III B 3, where we investigate effects of heterogeneity, the initial synaptic weight is assumed to be g 0 common for all the synapses. We vary g 0 as a control parameter.

Threshold for frequency synchrony
We compare how STDP affects the possibility of entrainment and formation of feedforward networks when a pacemaker is present and when absent. To this end, we fix a random network and a realization of ω i (1 ≤ i ≤ N). Without loss of generality, we assume  ω 1 ≥ ω 2 ≥ · · · ≥ ω N . For the network with a pacemaker, we make the fastest neuron a pacemaker. By definition, the rhythm of the pacemaker is not affected by those of the other neurons even though the pacemaker is postsynaptic to approximately k neurons. Using the bisection method, we determine the threshold value of g 0 above which all the neurons will synchronize in frequency.
The results without dynamical noise (i.e., σ = 0) are summarized in Table I. When the pacemaker is present from the beginning, STDP drastically reduces the threshold for entrainment [11]. After entrainment, all the neurons rotate at the inherent frequency of the pacemaker, that is, ω 1 = 8.60. When a pacemaker is initially absent, STDP reduces the threshold for frequency synchrony by 34%. Facilitation of frequency synchrony in the absence of the initial pacemaker is consistent with the results for the complete graph with N = 3 (Fig. 2). In this situation, the scenario to frequency synchrony is different between the presence and the absence of STDP. With STDP, the fastest oscillator eventually entrains the entire network when the initial synaptic weight is above the threshold, as in the case of the network with a prescribed pacemaker. Without STDP, the fastest oscillator does not entrain the other neurons. The realized mean frequency 8.08 is close to the mean inherent frequency of the 100 neurons. This suggests that frequency synchrony in this case is achieved by mutual interaction, rather than by one-way interaction underlying the entrainment by the fastest neuron. Therefore, in networks without predetermined pacemakers, STDP enables emergence of pacemakers and changes the collective dynamics drastically.

Network dynamics
For σ = 0, example rastergrams when there is initially no pacemaker and g 0 = 1.0, which is above the threshold value 0.72 (see Table I), are shown in Fig. 6. Figures 6(a) and 6(b) correspond to the initial and final stages of a simulation run under STDP, respectively; frequency synchrony appears as a result of STDP. Figure 6(c), which is an enlargement of Fig. 6(b), shows that the fastest neuron entrains the other neurons and that faster neurons tend to fire earlier in a cycle. Figure 7 shows the time course of the degree of synchrony r. Around t = 1.2 x 10 7 , r sharply drops, and all the neurons start to oscillate at the same frequency. The effective network defined by the surviving synapses in the final state is drawn in Fig. 8. The neurons are placed so that the horizontal position represents relative spike time in a cycle. With this ordering, the neurons form a feedforward network. In other words, after STDP, if a presynaptic neuron fires later than a postsynaptic neuron in a cycle, this synapse is not present.
Partial entrainment occurs when g 0 is slightly or moderately smaller than the threshold value 0.72. Circles and crosses in Fig. 9 represent the actual frequency after transient and the inherent frequency of the each neuron, respectively, when g 0 = 0.5. The neurons with the same actual frequency belong to the same cluster. Each cluster forms a feedforward network emanating from an emergent pacemaker. Figure 9 indicates that the neurons are divided into two clusters and one isolated neuron. Neuron 2 entrains 85 other neurons all of which are slower than neuron 2, neuron 6 entrains 12 slower neurons, and neuron 1 gets isolated. In this and further numerical simulations we performed, the root of a feedforward subnetwork is always occupied by the fastest neuron in the cluster.
Whether two neurons eventually belong to the same cluster is determined by where these neurons are located on the initial random network and by how close their inherent frequencies are. If g 0 is smaller than the value used for Fig. 9, two neurons have to be closer in ω i to stay connected after STDP. Then the number of clusters increases, and the number of neurons in a cluster decreases on average.

Robustness against dynamical noise and heterogeneity
To examine the robustness of the numerical results reported in Sec. III B 2, we perform additional numerical simulations with dynamical noise and random initial synaptic weights.
We draw initial g ji (j, i) ∈ E randomly and independently from the uniform density on In sum, emergence of entrainment via STDP survives some dynamical noise and heterogeneity in the initial synaptic weights. We have confirmed that, when the entrainment occurs, it is quickly established at around t = 10 6 , and the fastest oscillator is located at the root of the feedforward network, as in Fig. 8.

Network motifs
We investigated the evolution of three-neuron networks in Sec. III A because we expect that these results have something common with evolution of such subnetworks in large networks. The results in Sec. III A predict the following: • Bidirectional edges do not survive STDP, and feedforward networks of size three will be relatively abundant after STDP. Subnetworks abundant in a large network relative to the case of random networks with the same mean degree (or other order parameters) are called network motifs [22]. The hypothesis that feedforward networks are motifs in large neural networks is consistent with the observations in C. elegans neural networks [22].
• As a result of STDP, a neuron has at most one effective upstream neuron unless multiple upstream neurons are very close in frequency.
There are 13 connected network patterns of three nodes. How often each pattern appears in a network with N = 100, relative to the random network, can be quantified by the Z score [22]. The Z score is the normalized number of a pattern in the network, where normalization is given by the mean and the standard deviation of the count of the pattern based on independent samples of the randomized networks. A pattern with a large Z score is a motif of the network with N = 100. Figure 11 shows the Z score of each pattern before (circle) and after (square) STDP, calculated by m-finder [23]. We set σ = 0 (i.e., no dynamical noise) in this analysis. The error bar shows a range of one standard deviation based on ten simulation runs in each of which we draw a different initial random network and a different realization of ω i (1 ≤ i ≤ 100).
Before STDP, the neural network is a directed random graph, so that the Z score for each pattern is around zero, meaning that no pattern is overrepresented or underrepresented significantly. After STDP, the feedforward network whose emergence and survival were observed in Sec. III A (i.e., pattern 5 in Fig. 11) and patterns consistent with this (i.e., patterns 1 and 2) are overrepresented. These are motifs of our final networks. Pattern 4 is also a motif in spite of our negation in Sec. III A because the two upstream neurons in pattern 4 have the same actual frequency. They are generally different in inherent frequency but share a more upstream ancestor. As the example network in Fig. 8 shows, existence of multiple paths from a neuron to another due to branching and uniting of edges is compatible with STDP. The other network patterns are not significant or underrepresented. These results are further evidence that feedforward networks are formed by STDP in heterogeneous neural networks.

IV. DISCUSSION
We have shown using heterogeneous coupled phase oscillators that feedforward networks spontaneously emerge via STDP when the initial synaptic weights are above the threshold value. When this is the case, the pacemaker emerges at the root of the feedforward network and entrains the others to oscillate at the inherent frequency of the pacemaker. Although these results have been known for two-neuron networks [11,13], we have shown them for the cases of three and more neurons and quantified the phase transitions separating frequency synchrony and asynchrony. The route to frequency synchrony is distinct from a conventional route to frequency synchrony that occurs when mutual, but not one-way, coupling between oscillators is strong enough. Some results obtained in this work are unique to the networks without a prescribed pacemaker. First, the emergent pacemaker is the fastest oscillator neuron according to our extensive numerical simulations. Note that all the oscillators fire at this frequency in the entrained state, whereas they fire at the mean inherent frequency of the oscillators when the frequency synchrony is realized by strong mutual coupling in the absence of STDP. Second, when the initial coupling strength is subthreshold, the neurons are segregated into clusters of feedforward networks. Third, our numerical evidence suggests that entrainment under STDP occurs more easily when a prescribed pacemaker is absent than present.
In spite of a wealth of evidence that real neural circuits are full of recurrent connectivity [24], feedforward structure may be embedded in recurrent neural networks for reliable transmission of information [12,25]. Feedforward transmission of synchronous volleys in rather homogeneous neural networks as those used in this work serves as a basis of reproducible transmission of more complex spatiotemporal spike patterns in more heterogeneous networks. Such patterns may code external input or appear as a result of neural computation [7,12]. Feedforward structure is also a viable mechanism for traveling waves often found in the brain [26]. Although computational roles of feedforward network structure are not sufficiently identified, our results give a support to the biological relevance of feedforward networks. The formation of feedforward networks, which we have shown for oscillatory neurons, is consistent with numerical results for more realistic excitable neurons subject to STDP [27]. The neurons that directly receive external input may be more excited and fire at a higher rate compared to other parts of a neural circuit. Our results suggest that such a neuron or an ensemble of neurons is capable of recruiting other neurons into entrainment and creating feedforward structure.
We assumed the additive STDP with the nearest-neighbor rule in which the dependence of the amount of plasticity on the current synaptic weight and the effects of distant presynaptic and postsynaptic spike pairs, triplets, and so on, are neglected. Generally speaking, evolution of synaptic weights is affected by the implementation of the STDP rule [18]. However, we believe that our results are robust in the variation of the STDP rule as far as it respects the enhancement of causal relationships between presynaptic and postsynaptic pairs of neurons.
Our preliminary numerical data with excitable neuron models suggest that the results are similar between the multiplicative rule [18] and the additive rule [Kato and Ikeguchi (private communication)]. Recent reports claim the relevance of acausal spike pairs in the presence of synaptic delay [10,28]. This and other factors, such as different time scales of LTP and LTD [27], may let bidirectional synapses survive as observed in in vitro experiments [29].
Incorporating these factors is an important future problem.
We have ignored inhibitory neurons for two reasons. First, our main goal is to identify phase transitions regarding synchrony with a simple model. Second, specific rules of STDP are not established for inhibitory neurons, albeit some pioneering results [30]. Taking inhibitory neurons into account, preferably in the subthreshold regime, warrants for future work. a simulation run. Similarly, the green is added when g 23 survives, and the blue is added when all the neurons get disconnected. Yellow regions appear when both g 12 and g 23 survive, since the combination of red and green is yellow. In this case, it turns out that g 13 also survives. We verified that no other connectivity, such as survival of g 13 without survival of g 12 or g 23 , and survival of g 21 , g 31 , or g 32 , appears except at points near phase transitions and resonance. Near phase transitions, we exclude such exceptional runs from the statistics. In the resonance regions (e.g., ∆ω ≈ 2.7 and g 0 ≈ 0.4), the three neurons may remain connected. In this situation, however, synaptic weights   is an enlargement of (b). We set N = 100, g 0 = 1.0, and σ = 0. The neurons are aligned according to the order of the inherent frequency.