Oscillation, Conduction Delays, and Learning Cooperate to Establish Neural Competition in Recurrent Networks

Specific memory might be stored in a subnetwork consisting of a small population of neurons. To select neurons involved in memory formation, neural competition might be essential. In this paper, we show that excitable neurons are competitive and organize into two assemblies in a recurrent network with spike timing-dependent synaptic plasticity (STDP) and axonal conduction delays. Neural competition is established by the cooperation of spontaneously induced neural oscillation, axonal conduction delays, and STDP. We also suggest that the competition mechanism in this paper is one of the basic functions required to organize memory-storing subnetworks into fine-scale cortical networks.


Introduction
In mice experiments [1,2], a memory is recalled when neurons that are active during a learning process are activated with optogenetic stimulation. A specific memory is considered to be stored in a subnetwork consisting of a small population of neurons. For such memory formation, competition among neurons might be necessary to embed memory-storing subnetworks into neural circuits [3][4][5][6]. Furthermore, synaptic plasticity is thought to play a critical role in subnetwork organization [7,8].
In the last decade, studies in the field of neuroscience have revealed that synaptic modification depends on presynaptic and postsynaptic neuronal activities; this is called spike timingdependent synaptic plasticity (STDP) [9][10][11][12]. Experimental observations suggest that the precise timing of presynaptic and postsynaptic neuronal action potentials is a crucial factor in information processing and/or memory formation in the brain. Therefore, STDP is thought to be one of the mechanisms to encode information into patterns of their synaptic weights [13,14].
Based on experimental observations, several studies have proposed models of STDP window functions [15][16][17][18]. For example, Song et al. [15] modeled synaptic learning that is independent of the synaptic weight; this is called an additive model or a hard-bound model. In contrast, Rossum et al. Rossum et al. [16] and Rubin et al. [17] proposed a model that linearly depends on the synaptic weight; this is called a multiplicative model or a soft-bound model. The

Neural network model
We employ Izhikevich's simple neuron model as the basis of our neural network [31,32]. This model is not only computationally effective as the leaky integrate-and-fire model, but can also realize firing patterns as rich as those in the Hodgkin-Huxley model [33]. Dynamics of the jth (j = 1, 2, . . ., N) neuron is described by the following two-dimensional ordinary differential equations: _ v j ¼ 0:04v 2 j þ 5v j þ 140 À u j þ I j ðtÞ; ð1Þ where v j is the membrane potential and u j is the recovery variable of the jth neuron. The membrane potential and the recovery variable of the neuron model is reset to c j mV and u j + d j when v j reached 30mV. The variable I j (t) represents inputs to the jth neuron at time t. The inputs are the summation of external inputs (I ext j ðtÞ) and synaptic inputs (I syn j ðtÞ). For the sake of simplicity, we model the synaptic inputs with the delta function δ(Á), where w ij is the synaptic weight from the ith neuron to the jth neuron, and t i, k is the arrival time of the kth (k = 1, 2, . . ., n j, i where n j, i represents the number of spikes of the ith presynaptic terminal of the jth neuron) spike at the ith presynaptic terminal. Our neural network consists of N (= 1,000) neurons including both excitatory and inhibitory neurons. The ratio of the excitatory neurons to the inhibitory neurons is 4: 1 [34]. In this paper, we use regular-spiking excitatory neurons and fast-spiking inhibitory neurons. The parameters for the excitatory neurons are set as a j = 0.02, b j = 0.2, c j = − 65, and d j = 8 and those for the inhibitory neurons are set as a j = 0.1, b j = 0.2, c j = − 65, and d j = 2 [31,32]. These neurons are randomly connected.
Because, in many experimental studies, the connection probability has been estimated between 0.1 and 0.3 [35][36][37][38][39], we choose 0.1 connection probability. In the neural network, no connections exist between any pairs of the inhibitory neurons. Furthermore, no neurons are self-connected. Excitatory connections have conduction delays of 1 to 10 ms, with a uniform distribution [40]. A time of 1 ms is required to transmit spikes on all inhibitory connections [40]. It has been shown that dendritic delays tend to strengthen self-feedback, whereas axonal delays weaken it [41,42]. We assume that the conduction delays are only axonal. Each neuron in the network receives an independent and uncorrelated Poisson spike train with the fixed firing rate of f spk/s during our simulation through a non-plastic excitatory feedforward connection. The spike train is statistically identical for both excitatory and inhibitory neurons. In the simulations, we test f = 1, 10, and 40 spk/s. The reason for the usage of the Poisson spike train is based on the observation that in vivo neuronal behaviors in cortical areas are highly irregular [29]. The amplitude of each spike in the spike sequence is set to 20 mV. In other words, I ext j ðtÞ ¼ 20 in Eq (1), which corresponds to a suprathreshold input when a neuron is in the resting state. A spike train for a neuron is statistically equivalent to spike trains for the other neurons. All excitatory synaptic weights are initially set to 6 mV, except for the simulation in the section of Independence of initial distribution of plastic synapses on neural competition. whereas all inhibitory synaptic weights are set to − 5 mV. Research has shown the significance of synaptic types [12], and therefore STDP is applied only to the excitatory synapses between excitatory neurons in the network.
STDP is a type of Hebbian synaptic plasticity that has attracted considerable attention [9,11,12]. In this synaptic plasticity, if a postsynaptic action potential follows a presynaptic action potential within tens of milliseconds, the synaptic weight between them is strengthened; this is long-term potentiation (LTP). On the other hand, a synapse is depressed if a presynaptic action potential follows a postsynaptic action potential, which is long-term depression (LTD). In this paper, we adopt the additive STDP rule proposed by Song et al. [15]. Its window function is expressed in terms of the exponential functions as follows: ( where Δt = t j − t i is the relative spike timing between a presynaptic terminal and a postsynaptic neuron. Hard bounds is assumed for plastic synapses. Therefore, the plastic synapses are constrained in the range of [0, w max ], where w max is set to 10 (except for the simulation in the section of Influence of neural network parameters on neural competition). The variable λ(= 0.1 mV) is the learning rate [40]. The variable α is the degree of asymmetry between LTD and LTP. This parameter is typically set to 1.2, but is varied in the simulation in the section of Influence of neural network parameters on neural competition. We use the same time constant τ (= 20 ms) for both the LTP and the LTD [12]. Synaptic derivatives are changed at individual firing events, and actual synaptic weights are updated once a second. In all numerical simulations, spike interactions in the STDP rule are limited to nearest-neighbor pairs, except for the simulation in the section of Influence of spike interactions in STDP on neural competition. [42]. Potentiation and depression, which are independent of firing events, are also included in the synaptic modifications as in Izhikevich [40].

Strength correlations
To quantify network structures, the degree of individual nodes is usually measured. If connections in networks are directional, we can take into account two types of degrees: indegree and outdegree. The indegree and the outdegree of the jth node can be defined by the total numbers of incoming (afferent) and outgoing (efferent) connections, respectively, and are expressed as where H 0 (x) is the Heaviside step function in which A high indegree implies that a neuron is affected by many other neurons, whereas a high outdegree implies that a neuron affects many other neurons through synaptic connections. If indegrees and outdegrees of neurons are biased, the bias is visualized in a joint degree distribution matrix (JDDM). The imbalance of indegrees and/or outdegrees of neurons in a network appears in distances from the main diagonal of the matrix. Degree distributions in real networks often have the scale-invariant or scale-free property [43].
In the case of neural networks, synaptic connections do not have only directions, but also weights. The quantification of such weighted directional networks needs a natural extension of the degrees defined in Eqs (5) and (6) [44]. These are called instrength and outstrength. The instrength and the outstrength of the jth neuron are defined by the sum of the normalized synaptic weights of afferent and efferent connections, respectively: An instrength indicates how much a neuron is affected by other neurons, whereas an outstrength indicates how much a neuron influences the other neurons. The imbalance between instrengths and outstrengths of neurons in a network is visualized in a joint strength distribution matrix (JSDM).
In the following of this study, we quantify self-organized neural network structures through STDP using the instrength and the outstrength defined by Eqs (7) and (8). Only excitatory connections between excitatory neurons in our network are plastic; therefore, we focus on a subnetwork consisting of the excitatory neurons for the network structure quantification.
To evaluate connectivity among neurons in our neural network, we introduce measures that we call instrength-and outstrength-correlation coefficient. The correlations are an extension of the degree correlation [45]. The degree correlation coefficient is usually computed from the total degree of the jth neuron: To calculate the degree correlation coefficient, remaining degrees are used. The remaining degree is the number that is one less than the total degree. However, in our calculation, instrengths and outstrengths are directly used. We define the instrength-and the outstrength-correlation coefficient as follows: where M = P i;j (See also the section of Connectivity of winner and loser neurons). The subscripts t and h denote the tail and the head of a connection and the superscript a indicates in or out.

Identification of winner or loser neurons
To identify the composition of a winner and a loser group, winner and loser neurons are defined based on their instrength. This is because, as seen in Fig 1, the instrengths tended to spread more widely than the outstrengths. As such, the threshold of the winner and the loser neurons is able to be easily determined. It should be noted that even if we have used the outstrength as a standard, our results will have been consistent in principle. The method to identify the winner groups and the loser groups is as follows. First, neurons are organized in descending order of their instrengths, and a neuron that satisfies the condition of s in < s out is identified. This neuron is set as the threshold and neurons before it are treated as winner neurons, and the remaining neurons are treated as loser neurons.

Definition of phase
To characterize activities of neural networks or populations of presynaptic terminals, we also obtain phases from oscillatory firing rates. These firing rates are low-pass filtered to define the phases, and the cut-off frequency of the filter is determined to be 35 Hz based on the result of our spectral analysis in the section of Neural competition is organized in neural oscillation. After filtering, we define the phases [46] as where t k and t k + 1 correspond to any pairs of neighboring negative peaks of the oscillatory firing rate. The subtraction of π is to arrange the positive peak between t k and t k + 1 at ϕ(t) = 0. We obtain the phases from the firing rates of the entire network or the presynaptic terminals.

Estimation of firing rates
The mean firing rate of N neurons is estimated by the following equation: where f j (t) is the firing rate of the jth neuron at time t. The firing rate of each neuron is given by where T ( = 10 ms) is the width of a temporal window, and ρ j (t) is the spike train of the jth neuron: where n j is the number of spikes of the jth neuron and t j, k denotes the time of the kth spike of the jth neuron. If the neuron emits a spike ρ j (t) = 1 and otherwise ρ j (t) = 0. In the case of presynaptic terminals, we only consider the conduction delays. That is, we use the timings of the presynaptic-terminal firings instead of the timings of the somatic firings and then conduct the same estimation.

Kendall's correlation coefficient
Kendall's correlation coefficient τ K is a non-parametric value to quantify correlation of a paired data set. Let us define a pair of data as x i and y i (i = 1, 2, . . ., m). We describe their ranks as X i and Y i and consider pairs of rank data (X i , Y i ). The data pairs are sorted in ascending order of Using the values P i and Q i (i = 1, 2, . . ., m − 1), Kendall's correlation coefficient is computed as where − 1 τ K 1.

Watson's U 2 -test
Watson's U 2 -test is a non-parametric test for phase data. This test identifies significant differences of the mean value and/or the variance of phase distributions. Here we describe ϕ x (i) (i = 1, 2, . . .m 1 ) and ϕ y (j) (j = 1, 2, . . ., m 2 ) as two samples where both are sorted in ascending order. Therefore, the indices of i and j represent the ranks of phase data in ϕ x (i) and ϕ y (j), respectively. The total amount of data is M = m 1 + m 2 . Next two variables X k and Y k (k = 1, 2, . . ., M) are prepared. For all the data of ϕ x (i) and ϕ y (j), the following process is repeated: In the table of significant values, the significance of the mean value and/or the variance in two-phase distributions is evaluated.

Equilibrium states of distributions of plastic synapses
First, we show synaptic distributions to check the temporal behavior of the synapses (Fig 1A-1C). In agreement with the previous numerical and theoretical studies with the Fokker-Plank theory [15,17,18,47], it is observed in our simulations that plastic synapses in the network bimodally distribute after the long-time simulations (t = 3,600 s). Due to the ability of the STDP to prevent a firing rate in networks from drastically increasing, many plastic synapses go to the lower bound if the firing rate in the networks is high [15,17,18]. Also after t = 3,600 s, individual synapses continually change due to firing events. Nevertheless, the influence of these changes are small and trivial, and the form of the synaptic distributions is almost invariant (S1 Fig). Then, we regard the networks after 3,600 s as being enough stable to quantify the network organization using synaptic weights. The stability of the organized network is further discussed in the following section.

STDP induces neural competition: emergence of winner and loser neurons
The joint strength distribution matrix (JSDM) represents the imbalance between the instrengths (s in ) and the outstrengths (s out ) observed in our neural network (Fig 1D-1F). The JSDM exhibits a two-dimensional Gaussian distribution before the STDP learning (results are not shown) because excitatory synapses between excitatory neurons are randomly connected and the weights of the synapses are homogeneous. After the STDP learning, the outstrengths are widely distributed, while the instrength distribution is narrow when the firing rates of the external inputs are 1 spk/s (Fig 1D).
For external inputs of 10 spk/s, the bias of the instrengths is magnified (Fig 1E). Moreover, the neurons in the neural network compete and two peaks emerge, indicating the existence of two assemblies. The neurons in one assembly have high instrength but low outstrength, while the neurons in the other assembly exhibit an opposite trend. The instrengths of the neurons in these two assemblies display clear differences. The outstrengths are also widely distributed but are narrower than those with 1 spk/s external input (Fig 1D and 1E).
In the case of external inputs with a higher firing rate of 40 spk/s, almost all excitatory neurons have similar instrengths and outstrengths, however, a few neurons achieve high instrength ( Fig 1F). In comparison to the outstrengths, the instrengths form a relatively wide distribution in the space of s in -s out .
In all the cases, synaptic competition is observed, but the ratio of synapses reaching the upper-bound and the lower-bound depends on the mean firing rate of the external inputs ( Fig  1A-1C). Additionally, the degree of neuronal competition changes depending on the mean firing rates (Fig 1D-1F). According to these results, the neuronal competition is related to the ratio of depressed synapses to potentiated synapses.

Stability of neural competition through STDP
To evaluate the stability of the neural competition, we count the number of neurons that move between the winner and the loser assembly from t to t + 1 s (Fig 2). In both cases, i.e. movement from the winner to the loser, or from the loser to the winner assemblies, a maximum of three neurons moved, corresponding to 0.38% of the total excitatory neurons. This change is negligible on a whole network. We, therefore, regard the neural competition as stable.

Connectivity of winner and loser neurons
To analyze how the competitive neurons are connected in the network, we characterize the neural network with instrength-and outstrength-correlation coefficient, r in and r out . These coefficients quantify the similarity of neurons at the ends of connections in networks (See also Fig 3A). When two neurons at the ends of connections in a network tend to have similar instrengths or outstrengths, the coefficients are positive. When instrengths or outstrengths are dissimilar, these coefficients are negative.
The time courses of r in and r out are shown in Fig 3B and 3C. At t = 0 s, both r in and r out are zero in all cases because the neurons in the neural network are randomly connected through synapses under the initial condition. Evidently, r in and r out only significantly decrease from the  initial condition at 10 spk/s. They then converge at approximately − 0.1 and − 0.12 (P < 0.001, t-test), respectively. In the other cases, r in and r out do not reach this level of dissimilarity. Such significant differences of the coefficient values for 10 spk/s case come from the clear competition between the winner and the loser assembly (Fig 1E). We should note that it is impossible for the coefficients to be negative at the significant level if the neurons do not compete and if the neurons in the individual assemblies have the reverse trend of the instrengths and the outstrengths as Fig 1E. Taking into account the results in Fig 1, the neurons are competitive and the dissimilar neurons tend to locate at the ends of individual connections. In other words, connections between the winner and the loser assembly are more strengthened and internal connections in the individual assemblies are easily pruned off.
In addition to the dismilarity, we also analyze the small-world property of the neural network [48]. At the beginning of the simulation (t = 0 s of Fig 2 in [48]), the network has as a large clustering coefficient and a characteristic path length as regular networks. The effect of the STDP leads to a much smaller characteristic path length relatively to the regular networks, although the large clustering coefficient is maintained (t = 600 s of Fig 2 in [48]). This indicates that the small-world network emerges in the connectivity among the winner and the loser neurons. In this study, neural competition is primary interest and the analyses therefore focus on cases of 10 spk/s-external input.

Independence of the numbers of excitatory and inhibitory presynaptic neurons on neural competition
We have shown that neurons in the network are competitive when 10 spk/s-external inputs are given to each neuron. Since presynaptic and postsynaptic activities determine synaptic modifications, the factors determining postsynaptic activity might be the key to neural competition. When considering how to construct our network (See Materials and Methods), one of the differences of individual neurons is the number of presynaptic neurons because we have adopted statistically equivalent external inputs. Therefore, for individual neurons receiving 10 spk/s external inputs, we plot the number of excitatory presynaptic neurons against their instrengths after learning (t = 3,600 s) as shown in Fig 4A. They clearly do not have a linear correlation so their correlation is quantified by the Kendall correlation coefficient τ K -a non-parametric index (See details in Materials and Methods). However, even this coefficient cannot identify a correlation (τ K = 0.18).
In the same way, the number of inhibitory presynaptic connections is plotted against the instrengths after STDP (t = 3,600 s) in Fig 4B. Since inhibitory synaptic weights are negative, they have a negative correlation, but the coefficient is small and negative (τ K = − 0.3). Therefore, we conclude that the number of excitatory and inhibitory synapses on individual neurons has a minimal effect on whether neurons obtain many or few strong synapses.

Independence of initial distribution of plastic synapses on neural competition
For further analysis of the independence of the initial network architecture, we also change a given distribution of excitatory synaptic weights in the initial condition and observed the JSDMs after STDP (t = 3,600 s) for the 10 spk/s external input samples in S1 File. The variety of the initial synaptic distribution does not affect the JSDM after learning, and STDP induces neuronal competition. The three JSDMs are in perfect agreement with Fig 1B, which indicates that the structures are robust for certain synaptic weights under the initial condition. Therefore, synaptic weights before learning has little influence on neuronal competition.

Relation between axonal conduction delay and synaptic modifications
Another conceivable difference in the statistics of presynaptic connections for individual neurons is the distribution of axonal conduction delays. As such, we analyze the ratio of each conduction delay to the total number of presynaptic connections for each neuron. Fig 5A shows the mean values of this ratio in the winner and the loser assembly. The smaller the conduction delays are, the higher their ratio is in the winner neurons. In contrast, the ratio of larger conduction delays is higher in the loser neurons. From this result, the high ratio of shorter-delay connections seems to be advantageous for neurons to become winners.
We also analyze the relation between the conduction delays and the mean synaptic weights after learning (t = 3,600 s) of both winner and loser neurons (Fig 5B). For both winner and loser neurons, the smallest delay has the largest mean weight, and the mean weight decreases as the conduction delay increases. This result indicates that smaller conduction delays lead to stronger synapses. However, during the first 5 ms of the conduction delays, the mean synaptic weights of the winner neurons are twice as strong as those of the loser neurons. This explains why winner neurons can get a high instrength but loser neurons cannot.

Neural competition is organized in neural oscillation
We have shown that winner neurons have a higher ratio of presynaptic connections with smaller conduction delays that are capable of more potentiation. However, it is still unclear why the small delay connections of the winner neurons are more strengthened than those of the loser neurons. To unveil origins of strong potentiation of the small delay synaptic connections of the winner neurons, we observe the network activity in the STDP network because it has a strong effect on synaptic modifications.
To show changes of the network activity during the learning process, we pick up the activity for 0-0.5 s, 3-3.5 s, and 5-5.5 s (Fig 6A-6C). The neural network exhibits oscillatory behaviors with frequency variations over time. At the beginning of the learning process, the network exhibits slow oscillation (Fig 6A), which speed up over time (Fig 6B and 6C). Additionally, the mean excitatory synaptic weight gradually decreases due to STDP and normalizes at around 2 ( Fig 6D).
To gain further understanding, we conduct a spectrum analysis on these oscillations ( Fig  6E). We notice that the time course of the power spectrum of the excitatory population looks very similar to that of the inhibitory population. This is explained as follows: By receiving the external inputs, excitatory and inhibitory neurons begin to fire and increase their firing rates. Triggered by the increase of the firing rate of the excitatory neurons, the inhibitory neurons are strongly activated. Indeed, the local maxima between the oscillations of excitatory and inhibitory neurons exhibit a small gap (Black arrows in Fig 6C bottom panel). Excessive firing of inhibitory neurons inactivates excitatory neurons by large amounts of negative feedback. This inhibition also leads to the silencing of the inhibitory neurons because of diminished excitatory inputs. The excitatory and the inhibitory neurons, however, are excited by the external inputs and begin to fire again. This cycle is repeated in the neural network, and therefore, this cycle results in the stable oscillation. This looks like the phenomenon known as the pyramidal-interneuronal gamma [49][50][51][52].
For comparison, we also observe the network behavior for the 1 spk/s-external input (Fig 7), where neural competition did not emerge (See also Fig 1D). Analogously to the 10 spk/s case, the network exhibits the slow oscillation at the early stages of learning (Fig 7A and 7B). As learning progress, the oscillation vanishes and the neuronal firing rate diminishes ( Fig 7C). As seen in Fig 7E, no apparent peaks exist in the power spectra. Even though the settings are different than in the simulation using 10 spk/s-external input, the mean synaptic weight declines and normalizes at same point as in the case of 10 spk/s-external input (Fig 7D). A case comparison suggest that the neural oscillation supports enhancing synaptic potentiation more on the winner neurons than on the loser neurons.

Mechanisms of neural competition
To understand the competition mechanisms, we observed presynaptic and postsynaptic activities of a winner and a loser neuron randomly selected from individual populations (Fig 8A). The mean firing rates of the presynaptic terminals of the sampled winner and loser neurons oscillate in a similar manner as the firing rates of the excitatory neurons. However, the oscillation of the presynaptic terminals is slightly delayed from the oscillation of the excitatory neurons because of the axonal conduction delays (Fig 8A asterisks).
To quantify the delays, we evaluated gaps in the phases of the presynaptic oscillations from the excitatory oscillation as illustrated in Fig 8B. The mean probabilities of the phase gaps are plotted in Fig 8C. For both types of neurons, the positive gaps are collected, indicating that the presynaptic oscillations are delayed by the excitatory oscillation. The probability of phases gaps [0, π/4] in the winner neurons (red bars) is always higher than in the losers (blue bars). In contrast, the probabilities of the loser neurons exceeded those of the winner neurons in the larger phase range. This result is in good agreement with the previous result (Fig 5A) because the winner neurons have a higher ratio of smaller conduction delays. The distributions of the phase lags of the winners and the losers are significantly different (P < 0.001, Watson's U 2 -test). For the details of Watson's U 2 -test, see Materials and Methods. To understand the potentiation and the depression processes in the neural oscillation, it is also necessary to observe the behaviors of the winner and the loser neurons. An example of the time trace of the membrane potentials of the sampled winner and the sampled loser neuron is depicted in Fig 8A (Lower panel). These neurons are the postsynaptic neurons of the synapses shown in Fig 8A (Upper panel). The winner neuron tends to fire just after the local maxima of its presynaptic firing rate (Black arrows for the red line in Fig 8A), while the loser neuron seems to emit spikes before the peaks of its presynaptic firing rate (Black arrows for the blue line in Fig 8A) as depicted in Fig 8D. Here, we characterized spikes of the winner and the loser neurons by using the phases of the presynaptic oscillations. The probability distributions of the The lower panel is the time courses of the membrane potential of (red) a winner and (blue) a loser neuron. The upper panel is the average firing rates of (black) all excitatory neurons, (red) presynaptic terminals of the winner neuron, and (blue) presynaptic terminals of the loser neuron. These average firing rages are estimated by Eq (11). For the black line, timings of somatic firings are used. The red and the blue line are estimated with firing timings of presynaptic terminals, at which delay lengths are added to somatic firing timings of presynaptic neurons of the winner and the loser neuron. The winner and the loser neuron are randomly picked up from the identified groups at t = 3,600 s. (B) Schematic of a phase difference of the presynaptic oscillations from the global oscillation. The black, red and blue lines represent the mean firing rate of all excitatory neurons, the presynaptic terminals of the winner neuron, and the presynaptic terminals of the loser neuron, respectively. (C) Phase distributions of the local maxima in presynaptic terminal firing rate for the winner (red) and the loser (blue) neurons. The local maxima for the presynaptic firing rates are characterized by estimated phases relative to phases of a global firing rate of the excitatory neurons with Eq (10). The data for 45 s (from 5 s to 50 s) is used for the estimation. Results are however not different when using data after t = 50 s. Watson's U 2 -test is used to test for significant differences in the distribution pattern between the winner and the loser neurons. This is the non-parametric test for phase data used to indicate any significant differences in the mean value or the variance (P < 0.001). (D) Schematic of the phase difference between the mean firing rate of presynaptic terminals and a postsynaptic firing in a winner neuron (left) and a loser neuron (right). The dashed and the solid lines represent the mean firing rate of presynaptic terminals and the postsynaptic potential, respectively. (E) Same as (C), but the local maxima of the presynaptic firing rates are replaced by the spikes of the winner (red) and the loser (blue) neuron firing rates (P < 0.001).
doi:10.1371/journal.pone.0146044.g008 phases are shown in Fig 8E. In both distributions, the largest peaks locate around π/2. In the range of [− π, 0], the probabilities of the loser neurons always exceed those of the winner neurons. The positions of the winners and the losers are reversed in [0, π/2]. This fact indicates that synapses on winner neurons are more potentiated and less depressed than those on loser neurons. The mean synaptic potentiation (or depression) between on winner neurons and on loser neurons is, on average, significantly different because there is significant difference in their probability distributions of phases (P < 0.001, Watson's U 2 -test). Strong synapses are then, developed more on winner neurons than loser neurons.

Influence of neural network parameters on neural competition
In the previous sections, we have showed that the cooperation of spontaneously induced oscillatory behaviors, axonal conduction delays, and learning interacted with each other and resulted in the neuronal competition. The behaviors of neurons and synapses must however be influenced by certain parameters in the network, which in turn affect neural competition. As such, we investigated the influence of various parameters on the neural competition.
The first two parameters are α determining the balance between LTD and LTP, and the inhibitory synaptic weight because they have strong impacts on the synaptic distributions [15,16,47,[49][50][51][52]. To evaluate organized neural network structures in these two parameter spaces, we used the instrength-and the outstrength-correlation coefficient. As shown in the results in Fig 3, when neural competition occurs, these coefficients become negative.
As expected, these parameters drastically changed the network connectivity (Fig 9). The network only established neural competition at 10 and 40 spk/s. In addition, the competition is observed for α 1.2. The parameter is in the suitable range of the experimentally observed STDP window function [12,53]. The influence of the inhibition level on the competition is not observed in our result.
The next parameter is the upper bound of plastic synapses. Any other parameters are the same as Fig 1B. In spite of the greater or the smaller maximum synaptic weight, neural competition does not emerge (r in = 0.002 and r out = 0.006 in Fig 10A, r in = 0.13 and r out = 0.05 in Fig  10B, r in = 0.3 and r out = 0.2 in Fig 10C). It can be considered that spike effects dramatically affects neural competition [54,55].

Influence of spike interactions in STDP on neural competition
In addition to certain parameters, STDP learning similarly influences neural competition in the network. In this section, we analyzed the influence of the spike interactions of STDP on neural competition. S2 File shows the JSDMs, when we simulated the neural network with the all-to-all interactions of STDP. All parameters of the network are the same as Fig 1. The neural competition is quantitatively and qualitatively identical to those in the case of the nearestneighbor interactions. These results are expected because the all-to-all spike interactions only increase the frequency of the synaptic potentiation and depression and do not affect the competition mechanisms shown in the previous sections. Hence, the spike interactions in STDP are not a key factor for neural competition.

Discussion
In this paper, we showed that excitable neurons in a recurrent network spontaneously competed and organized into two assemblies: winners and losers. Our analyses revealed that spontaneously induced neural oscillation, axonal conduction delays, and learning cooperated to establish neural competition.
In our numerical simulations, on average, STDP decreased synaptic weights in the network. At the same time, a certain level of neural activation spontaneously induced neural oscillation due to the existence of inhibitory neurons in the network. Around the local maxima of the oscillation, instantaneous firing rates of neurons in the network were high. If a postsynaptic neuron received inputs through many small delay connections, the inputs from the other neurons immediately arrived at the neuron. The spikes arrived before the next local oscillation maximum. Then, the presynaptic inputs could effectively induce firings of the postsynaptic neuron. This leads to potentiation of many synapses. Neurons that had high ratio of small axonal conduction delays were able to be winners. In contrast, the spike arrival at presynaptic terminals is delayed in neurons with many large delay connections. Therefore, many synapses failed to be potentiated. In addition, the LTD window is dominant in STDP, which already depressed the synapses in the network. For these reasons, the neurons with many large delay connections became the losers.
Because a specific memory is thought to be stored in a subnetwork consisting of a small population of neural circuits [1,2], the neural competition shown in this paper is an important property for memory formation [3][4][5][6]. The competition mechanism in this paper is applicable to real neural circuits because all key factors of competition are observable in neural circuits. Because the length of a axonal conduction delay is proportional to the distance between a pair of neurons, our result indicates that synaptic connections of closer neurons were strengthened by neural oscillations. In other words, the neurons organized locally dense but globally sparse neural circuits. This might be related to the distance-dependent high-order correlations of neuronal activities [7,8]. It is implied that highly nonrandom local connectivity is organized by the distance-dependent high-order correlations [36,37].
In our analyses, the neural competition is accomplished when certain parameters of numerical simulations were set within physiologically reasonable limits. For example, the firing rate of external inputs, the slight dominance of the depression in the STDP window, and the upper bound of plastic synapses. The parameters tested in this paper might be unique in different areas of the brain [9,56]. As shown in this paper, it is suggested that local circuits organize very differently in different brain areas.
Iglesias et al. conducted similar simulations and analyses using the leaky-integrate and fire units and the additive STDP, in which spike interaction is the nearest neighbor, without axonal conduction delays [23]. They extended their simulations and analyses by introducing neuron death with larger size of neural networks [57,58]. In their simulations of Ref. [23], they assumed the application of inputs to a fraction of the population, that shaped a bar column, that dynamically moved, in 2D lattice. In Fig 4 of Ref. [23], Eqs (5) and (6) were used in their analyses. In the result, the neurons exhibited clear competition in the excitatory neural population. The indgrees in two neural groups were much different but their outdegrees were the same level. However, our analyses with the degrees show a reverse trend. The main difference of our simulation settings from Ref. [23] is the way of the application of inputs and the existence of conduction delays. In our result, we do not observe large difference of indegrees between the two groups (S2 Fig). Rather, there is the gap of outdegrees between them. This difference might come from the effect of the way of the external stimulation and the axonal conduction delays, especially the latter one is important for the neural competition shown in the current study because of the competition mechanisms (See also the section of Mechanisms of neural competition). Then, it is considered that the competition in the STDP neural networks in the current study originates in the different mechanisms from ones of Ref. [23]. So, it is suggested that there are some possible mechanisms of the neural competition. Accordingly, the neural competition in STDP recurrent networks should be further analyzed.
Supporting Information S1 File. Influence of initial synaptic distributions on the organization of networks. Initial excitatory synaptic weights have a uniform distribution in the range of ( Figure A) [5,7], ( Figure B) [4,8], and ( Figure C) [3,9]. The panels show JSDMs. The instrength-and outstrength-correlation coefficients for each case are ( Figure A)