Delay-Induced Multistability and Loop Formation in Neuronal Networks with Spike-Timing-Dependent Plasticity

Spike-timing-dependent plasticity (STDP) adjusts synaptic strengths according to the precise timing of pre- and postsynaptic spike pairs. Theoretical and computational studies have revealed that STDP may contribute to the emergence of a variety of structural and dynamical states in plastic neuronal populations. In this manuscript, we show that by incorporating dendritic and axonal propagation delays in recurrent networks of oscillatory neurons, the asymptotic connectivity displays multistability, where different structures emerge depending on the initial distribution of the synaptic strengths. In particular, we show that the standard deviation of the initial distribution of synaptic weights, besides its mean, determines the main properties of the emergent structural connectivity such as the mean final synaptic weight, the number of two-neuron loops and the symmetry of the final structure. We also show that the firing rates of the neurons affect the evolution of the network, and a more symmetric configuration of the synapses emerges at higher firing rates. We justify the network results based on a two-neuron framework and show how the results translate to large recurrent networks.

Dendritic and axonal propagation delays within and between brain areas may assume different values and need not be identical. Dendritic delays are typically smaller than axonal delays 20 . However, dendritic delays may range from sub-millisecond to a few milliseconds 39,40 , e.g. Agmon-Snir and Segev 39 demonstrated that the dendritic delay in octopus cells is about 0.5 ms. Axonal propagation delays may range from 0.3 ms in thalamo-cortical connections 41 and 1 ms in cortico-tectal connections 42 to 20 ms in cortico-cortical connections 43 , and even more, up to 40 ms in cortico-thalamic circuits 42 .
The presence of dendritic and axonal propagation delays can regulate the emergent structures of STDP-driven plastic neuronal populations 20,32 . Propagation delays affect the function of STDP in two different ways. First, the total propagation delay τ = τ d + τ a , i.e., the sum of the dendritic τ d and axonal τ a propagation delays, determine the synchronization tendency of the coupled neurons and the time difference of spiking of the neurons [44][45][46][47][48][49][50] . Second, since the effect of pre-and postsynaptic neurons arrive at the synapse after axonal and dendritic delays, respectively, the difference of the two propagation delays ξ = τ d − τ a , changes the relative timing at the synapse. We recently showed 32 that by incorporating propagation delays in a plastic neuronal network model, bidirectional connections can be retained and potentiated despite the previously reported simulation results with conventional STDP in the absence of delays 33,34,38 .
In this study, we focus on the multistability of the structural dynamics of the network with respect to the distribution of the initial synaptic weights in coupled networks with delays. Multistability was previously reported with respect to the mean initial synaptic weight: The final mean synaptic strength mainly depends on the initial average weight 34,38 . Here we show that within the present framework not only the mean, but also the standard deviation of the initial connections crucially determines the final mean strength of the synapses. Moreover, the presence of bidirectional connections and the symmetry of the final structure is determined by the initial discrepancy of the synaptic weights. We show that the multistability of the network can be explained in a modular way, by the presence of different attractors for a two-neuron motif in the two-dimensional space of the initial connection strengths. We also show that the domain of attraction of each final structural state depends on the firing frequency of the neurons: At high firing rates the attraction domain for symmetric connections grows in expense of shrinking that of the unidirectional connections. This leads to a more symmetric final structure when the neurons are firing at higher rates. Intriguingly, all of these nontrivial phenomena are only seen when the propagation delays are incorporated in the formulation of the model: Without propagation delays, any initial preparation ends up with unidirectional connections regardless of the firing rate of the neurons and the initial weights of the synapses.

Results
Theoretical Framework. We consider two neuronal oscillators described in a phase reduced model. In this approach, the evolution equation for the relative phase of the two neuronal oscillators can be written as follows (for the derivation see Methods, Eq. (6)): 12 21 12 21 where χ = ϕ 2 − ϕ 1 is the phase lag, π ν Ω = Δ 2 that Δν = ν 2 − ν 1 is the frequency mismatch, and ψ = 2πντ is the rescaled delay with τ = τ d + τ a being the sum of dendritic τ d and axonal τ a propagation delays. The values for dendritic and axonal propagation delays are chosen based on the experimental observations in cortical areas 42,43 . The quantity Γ reflects the relative difference of the synaptic strengths in two opposite directions, where g ij with i, j = 1, 2 (j ≠ i) is the synaptic strength of the synapse connecting presynaptic neuron j to postsynaptic neuron i.
The synapses connecting the two neurons are subjected to classic pair-based STDP rule where the synaptic strengths g ij = g ij (t) are modified based on the delayed time lag between pre-and postsynaptic activity at synaptic site, i.e., ξ Δ ′ = Δ + t t ij ij 32 where Δt ij = t i − t j is the time lag in the absence of propagation delays and ξ = τ d − τ a is the difference between dendritic and axonal propagation delays. The synaptic strengths are updated by an additive rule at each step g ij → g ij + Δg ij , according to the following learning function (see Methods): where A + (A − ) and τ + (τ − ) are the learning rate and the effective time window of synaptic potentiation (depression), respectively, and sgn(Δt) is the sign function. The synaptic strengths are bounded in the range (g min , g max ) = [0.05, 1] by the hard bound saturation constraint 4,51 : The synaptic strengths are set to g min (g max ) once they cross the lower (upper) limit of their allowed range. We consider a balanced STDP profile 4,19 with equal potentiation and depression contribution, i.e. A + = A − and τ + = τ − , in order to extract the pure influence of propagation delay and oscillation frequency on the synaptic dynamics. We assumed that the changes in synaptic strengths are slow in comparison to the fast timescale of the system, therefore, the neurons remain phase-locked despite of eventual change in the synaptic strengths. Given the firing frequency ν and assuming that the frequency mismatch between the two oscillators is negligible, i.e. Ω  0, the fast timescale equation of the two-neuron motif (Eq. (1)) exhibits a stable fixed point which is given by Eq. (9), χ ψ = Γ − ⁎ tan ( tan ) 1 , and hence firing time difference of two neurons at the stable state will be Δt * = χ * /(2πν). These values could be inserted in Eq. (3) to give the evolution of the synaptic connections 32,52 . Bistability in the Two-Neuron Motif. In Fig. 1 we show the synaptic dynamics (illustrated by the vector field) overlaid on the the representation of the instantaneous time lag of the spiking neurons (shown by colors), for the different values of the initial synaptic strengths (g 21 (0), g 12 (0)). It can be seen that the basin of attraction for the three different final connectivity patterns (bidirectional, unidirectional and decoupled) is determined by the firing frequency ν and the difference between dendritic and axonal delays ξ. The results shown in Fig. 1 denote the theoretical prediction of different stable structural states that coexist for different combinations of delay and firing frequency. The colors indicate the stable fixed point of the time lag Δt * given by Eq. (9), and the vector field shows the changes of the synaptic strengths from Eq. (3). Given the initial values of the synaptic strengths, the instantaneous color-coded fixed point of the time lag determines the synaptic changes depicted by the vector field, and the subsequent values of the synaptic strengths. The corresponding trajectories of the time evolution of the synaptic strengths resulted from numerical experiments with three different initial values are shown by yellow solid curves in Fig. 1G and I. The numerical solutions practically follow the vector field directions predicted by the analytical results. The time course of the simulated synaptic strengths and the firing time lag are presented in Fig. 2. Figure 1, left and right columns illustrate that according to different initial arrangements of the connections, the two-neuron motif will end up with one of the bistable coupling regimes: Bidirectional/unidirectional coupling ( Fig. 1, left column, dendritic delay is greater), or decoupled/unidirectional coupling (Fig. 1, right column, axonal delay is greater). The basin of attraction of the symmetric coupling regimes (i.e. bidirectional or decoupled final configuration) grows with an increase of the firing frequency (Fig. 1, left and right columns, top to bottom), making it more probable that in the high-frequency firing regime the neurons are either bidirectionally coupled (Fig. 1G) or remain decoupled (Fig. 1I).
We also examined the effect of a frequency mismatch Δν = ν 2 − ν 1 , i.e. the difference between the frequencies of the two oscillators introduced in Eq. (1), on the evolution of the synapses. Expectedly, the frequency mismatch shifts the border of the basins of attraction in a way that the neuron with the greater firing rate more likely wins the competition in generating a stronger outgoing synapse and suppressing the reverse one (see Fig. 3), which is in accordance with previous studies 21,53 . However, for small frequency mismatches (e.g. Δν = 0.02 Hz in Fig. 3A-C) bistability can still be present, and all three final configurations may be achieved depending on other parameters, i.e., the difference in propagation delays and firing frequencies.
The relation between dendritic and axonal delays is crucial to the emergence of bistability. When dendritic and axonal delays are equal, ξ = 0, which is equivalent to ignoring both in the STDP implementation, no bistability  In the presence of frequency mismatch the output synapse of the neuron with higher frequency is more likely to be potentiated, while the reverse synapse is depressed. emerges and the final connectivity pattern is unidirectional, regardless of the initial synaptic strengths (see Figs 1 and 3, middle columns). This is consistent with the results of Babadi and Abbott 23 , obtained for a balanced STDP profile. Note, in our model the bistability and the possibility of the emergence of symmetric connectivity patterns are observed for symmetric STDP just because the propagation delays are explicitly considered (also see Madadi Asl et al. 32 ). The dependence of the final stable coupling regimes on the initial values of the synaptic strengths can be explained based on Eq. (9), , that gives the asymptotic time lag between pre-and postsynaptic spikes. Therefore, the instantaneous fixed point of the time lag (Eq. (9)) is determined by the relative synaptic strength Γ, the firing frequency of the oscillation ν, and the total delay τ. In turn, the instantaneous value of Δt * determines the evolution of the synaptic strengths according to Eq. (3). More specifically, the values of the synaptic strengths determine the relative synaptic strength Γ = |(g 12 − g 21 )/(g 12 + g 21 )| from Eq. (1), where the interaction of the relative synaptic strength with the fixed point of the time lag Δ ⁎ t ultimately determines the final stable coupling regime. The right panels of Fig. 4, show that the color-coded relative synaptic strength Γ is close to zero around the diagonal line (light color) and maximal at the corners (dark color). As shown in Fig. 4 (left panels) the basin of attraction of the two-neuron motif changes when Δ ⁎ t (blue curve) crosses a certain value of Γ (denoted by arrow). The bistability is between bidirectional (red) and unidirectional (orange) connections

Recurrent Excitatory Networks.
In the two-neuron motif the relative synaptic strength Γ , given by Eq.
(1), determines the final configuration of the connections (see Fig. 4). Small Γ favors more symmetric final structure, and bidirectional connections can emerge if the axonal delays are small. Intriguingly, in large networks the initial distribution of the connection strengths is the key factor determining the emergence of bidirectional connections: In an initially homogeneous network with narrow distribution of the synaptic strengths, bidirectional connections (i.e. two-neuron loops) are more likely to survive in the network, while heterogeneity in the initial synaptic strengths drives the network to an asymmetric final structure with a small number of bidirectional connections. In contrast, so far, the evolution of the network and the emergent structure was attributed only to the mean of the initial synaptic strengths 33,34,38 . In order to study the generalization of our framework to the network level, we took a fixed value of the mean initial synaptic strength and varied the width of the distribution of the initial synaptic strengths by changing the standard deviation σ g of the distribution. We constructed an initially fully connected network of N = 200 excitatory neurons. Synaptic strengths are modified based on the pair-based additive STDP rule given by Eq. (2). Based on the analysis for the case of a two-neuron motif, the dynamics of the network for a given relative delay (difference of axonal and dendritic delays) and firing frequency can be predicted by the initial distribution of the synaptic strengths. Figure 5 shows the results for the recurrent excitatory network with relative delay ξ = 0.2 ms, as in the case of short-range synaptic projections, in the high-frequency firing regime. The dendritic delay τ d = 0.5 ms is greater than the axonal τ a = 0.3 ms. According to Fig. 1G, our analysis of the two-neuron motif predicts a multistability regarding the number of bidirectional and unidirectional connections in the emergent network. Figure 5A, left panel shows the time evolution of the mean synaptic strength, and right panel illustrates the initial and final distribution of the synaptic strengths where the network architecture is dominated by either bidirectional loops (blue distribution) or unidirectional connections (grey distribution). Figure 5B, left panel shows the synchronization degree of the neuronal oscillators and right panel represents the number of two-neuron loops. Figure 5C shows the final coupling matrices for pre-and postsynaptic neurons. For smaller values of the standard deviation (e.g. σ g = 0.05, 0.08 in Fig. 5C), based on our two-neuron motif analysis, the network is driven to potentiate bidirectional connections. However, for greater values of the standard deviation (e.g. σ g = 0.10, 0.15 in Fig. 5C), in the emergent structure unidirectional connections have more chance to survive.
For large axonal delays, as in the case of long-range synaptic projections, a similar bistability is observed between unidirectional connections and a disconnected network. Figure 6 shows the results for ξ = −0.5 ms. The axonal delay τ a = 1.0 ms is greater than the dendritic τ d = 0.5 ms. Figure 6A, left panel shows the evolution of the mean synaptic strength in the network, and right panel represents the initial and final distribution of the synaptic strengths where the network structure tends to be decoupled (blue distribution) or dominated by unidirectional connections (grey distribution). Figure 6B, left panel shows the synchronization degree of the neuronal oscillators and right panel represents the number of two-neuron loops. Figure 6C shows the final coupling matrices for preand postsynaptic neurons. In the case with small standard deviation of the synaptic strengths (e.g. σ g = 0.05, 0.08 in Fig. 6C), most of the connections tend to be depressed. In contrast, for wider distributions (e.g. σ = . . for g min stabilizes the network structure at a point where all connections attain the minimum value. Setting the minimum synaptic weight at zero, destabilizes the disconnected network, and one of the synapses builds up to achieve a unidirectional link. Our choice seems reasonable based on the experimental results on STDP 3,5 . The emergence of a stable decoupled state in the phase-locked mode can, e.g., be of importance for the development of stimulation methods for neurological disorders designed to induce long-lasting, sustained therapeutic effects by shifting a neuronal network from phase-locked attractors to more favorable, decoupled attractors 38,54,55 . By decreasing the firing frequency and increasing the disparity of the neuronal firing rate the symmetric (bidirectional or decoupled) connections are less likely to survive as revealed by the two-neuron motif analysis (Figs 1 and 3). We extracted the probability (P asym ) of appearance of unidirectional connections in the two-neuron motif, by taking the ratio of the area in the (g 21 , g 12 ) space which leads to unidirectional connections to the total area of the same parameter space. Figure 7A shows that this probability is a decreasing function of frequency. To investigate the impact of this result on the evolution of the network, we defined the network asymmetry index C net which quantifies the average asymmetry of the connections between pairs of neurons (see Methods). Expectedly, the network asymmetry index shows a similar decreasing trend which indicates that the number of unidirectional connections in the network decreases with increasing frequency (see Fig. 7A). As shown in Fig. 7B, by increasing the frequency of the oscillations, network asymmetry index decreases, which is consistent with a similar decrease in P asym in the two-neuron motif, indicating that at higher frequencies the probability of the appearance of symmetric connections (two-neuron loops or decoupled pairs) increases. Also, based on the results of Fig. 3, a sufficiently pronounced inhomogeneity in the firing rate of the neurons can prevent the formation of strong bidirectional connections. Figures 8 and 9  the distribution of the firing rates. Figures 8C and 9C show the final coupling matrices for pre-and postsynaptic neurons. As shown in Figs 8 and 9, by increasing the standard deviation σ ν of the distribution of the firing rates, the number of loops in the final structure decreases. This is also illustrated in Fig. 7C where the asymmetry index C net increases with increasing σ ν . The probability P asym in the two-neuron motif also follows the same ascending trend. Figure 7D shows that by increasing the inhomogeneity in the initial distribution of the synaptic strengths with standard deviation σ g , unidirectional connections with high values of the asymmetry index are more likely to emerge. This increasing trend was also generally predicted by the probability P asym in the two-neuron motif based on our theoretical framework.

Discussion
STDP is conventionally known as a mechanism which generically breaks the structural symmetry of neuronal networks and underlies the formation of feedforward networks 4,19,24,25,56 . This result in recurrent networks contradicts the experimental observations in cortical circuits where bidirectional connections are frequent [26][27][28][29][30] . Bidirectional connections can be retained through alternative versions of conventional STDP: Triplet-based STDP 57 and STDP with a shifted profile 23 . Triplet-based STDP, which has been designed to conform the dependence of the experimentally observed synaptic changes to the firing rates 58 , can promote bidirectional connections in the high-frequency regime 57 . Pair-based STDP with a rightward shifted profile can also lead to bidirectional connections 23 , if potentiation dominates depression. The imbalance of the STDP profile in favor of potentiation in short spike-pairing times, may result in the retention of bidirectional connections when the neurons are loosely phase-locked 23 . In a more realistic modeling, considering forward and backward propagation delays, we have previously shown 32 that the widely accepted consequences of conventional STDP, mentioned above, can be challenged. Here we showed that in this framework, the evolution of the bidirectional connections and the symmetry of the emerging structure shows multistability and the degree of heterogeneity in the initial setting of the network determines the final structure of the network. Experimentally observed pair-based STDP parameters have been shown to be unbalanced: A + > A − and τ + < τ − parameters, e.g. propagation delays, on the dynamics of the system. For example, as noted above, in a potentiation dominated profile of STDP, bidirectional connections can emerge without incorporation of propagation delays in the model 23,31,35 . In the present study we chose a balanced STDP profile to highlight the effects which can be directly attributed to propagation delays. Previous studies ascribed the final stable connection regime in neuronal networks with synaptic plasticity without propagation delays to the initial mean coupling 33,34 . Here, based on the analysis presented for the two-neuron motif, we focused on the multistability of the network structure due to the disparity of the initial synaptic strengths. Our results show that for every pair of neurons, the final configuration of the synaptic connections follows a bistable dynamics governed by the initial mismatch of the strength of the synaptic connections in both directions, where the initial asymmetry in the structure (difference in the synaptic strengths) specifies the final synaptic configuration. The emerging configurations depend on the difference between dendritic and axonal delays: For larger dendritic delays, the final structures are unidirectional and bidirectional, while for larger axonal delays the two achievable states are unidirectional or fully decoupled neurons. In particular, our results highlight the importance of propagation delays, which are usually neglected in modeling studies of STDP: Without taking into account delays, the bistability cannot emerge, and all initial settings end up with a unidirectional configuration.
Our findings demonstrate the crucial role of the presence and the range of dendritic and axonal propagation delays in regulating the emergent structures of synaptic connections in plastic excitatory neuronal networks. We show that connections with experimentally measured values of short-range delays 41,42 and gamma band firing rates 58 might be favorable for strong bidirectional loops or unidirectional connections, depending on the initial distribution of the synaptic strengths (see Fig. 5). In contrast, in the case of long-range connections, a loosely connected network or unidirectional synaptic structures might emerge (see Fig. 6). Therefore, the difference of dendritic and axonal propagation delays is crucial for the selection of the network's ultimately emerging coupling regimes. In this way, qualitatively distinct connectivity patterns may emerge. The emergence of multistable coupling regimes can be of interest since such systems can be used in order to construct memory building blocks. Bistability of the synaptic strengths, arising from the positive feedback nature of STDP, results in the emergence of different stable attractor states of the synaptic connections 59 . Switching between different stable attractors enables neuronal networks to maintain and retrieve memories. Bistability of the network dynamics with recurrent connections mediated by STDP has been addressed for stationary 13 or oscillatory inputs 60 . Bistable dynamics can be characterized by the presence of two extreme stable fixed points of the synaptic strengths at g min and g max , imposed by the hard bound saturation constraint.The choice of the saturation constraints can have significant effects on the final distribution of the synaptic strengths; however, it does not crucially affect the weight dynamics 61 . Therefore, the main results presented in this study, i.e., development of different connectivity patterns can also be qualitatively obtained when soft bounds are imposed on the synaptic strengths.
The dependence of the emergent structure to the firing rate of the neurons was experimentally observed 58 , showing that pairing of spikes with both positive and negative time differences leads to potentiation in the high firing rate regime, which implicitly indicates that bidirectional connections can be promoted in this regime. While previous attempts with pair-based STDP failed to reproduce these results, our study shows that this can be achieved by considering propagation delays without having to consider triplets of spikes in the formulation 32 .
Here, we showed that the domain of attraction for each achievable final configuration of the synapses depends on the firing frequency of the neurons. For low firing frequency, the two synapses in the pairs of bidirectional connections evolve in different directions and unidirectional connections are more likely to emerge and survive in the network, as suggested by several studies on STDP [18][19][20][21][22][23] . However, by increasing the firing frequency, the simultaneous potentiation or depression of the bidirectional connections becomes more probable, so that the measure of the attraction domain for a symmetric configuration in the two dimensional space of the initial synaptic strengths expands. Accordingly, in recurrent neuronal networks the width of the distribution of the synaptic strengths and the firing frequencies determine the emergent structure of the network and, in particular, the number of two-neuron loops in the final structure. Furthermore, our results show that the symmetric depression of reciprocal connections between neurons in the phase-locked state is possible when favorable combinations of initial synaptic strengths, range of firing rates, and propagation delays are met. We showed that by considering a non-zero lower hard bound g min , the depression of reciprocal connections can be stabilized in order to construct a loosely connected network, as predicted by our theoretical framework. The possibility of a symmetric stable depression of reciprocal connections in the phase-locked mode may contribute to a further development of brain stimulation techniques that induce an anti-kindling, i.e., an unlearning of abnormally up-regulated synaptic connectivity and, in turn, abnormal synchrony 38 . Coordinated reset (CR) stimulation 62 , a desynchronizing multisite stimulation technique was successfully tested in preclinical 63,64 and clinical [65][66][67] proof of concept studies. However, our approach presented here, may lead to further improvements of brain stimulation techniques causing sustained therapeutic effects.

Methods
Pair-Based STDP Model. The neuronal oscillators are subjected to classic pair-based STDP rule where the synaptic strengths are modified based on the learning window function introduced in Eq. (2). The evolution of the synaptic strengths by the pair-based STDP rule can be calculated by taking the average of Eq. (2) over a period of the spiking neuron and smoothing, Supplementary Fig. S1; for details on the interplay between the temporal scale of the parameters with STDP models see 68 ), which is given by: 21 12 where A + (A − ) and τ + (τ − ) are the rate and the effective time window of synaptic potentiation (depression), respectively. T is the period of spiking and sgn(Δt) is the sign function. In the entire manuscript the profile of STDP is balanced by setting the parameters to A + = A − = 0.005, and τ + = τ − = 20 ms. Synaptic strengths are confined in the interval (g min , g max ) = [0.05,1]. It should be noted that hard bounds are imposed on the allowed range of synaptic strengths: The synaptic strengths are set to g min (g max ) once they cross the lower (upper) limit of their allowed range.

Phase Reduced Model.
Considering that the rate of synaptic change is small, and therefore, the changes in the synaptic strengths are negligible on the fast timescale of the system, the reduced averaged phase model for weakly pulse-coupled neuronal oscillators characterized by intrinsic frequency ω i = 2πν i and infinitesimal phase sensitivity Z(ϕ) can be written as follows 69,70 : where the neuronal oscillators are coupled via delayed connections of strength g ij with total delay τ = τ d + τ a . ϕ i is the phase of the i-th oscillator, ψ = ω i τ is the rescaled delay, and N is the number of oscillators. The neurons fire every time their phase passes multiples of 2π. In the model, we ignore the synaptic processing time, but the results are not affected by this assumption. Eq. (4) can be written for two coupled type-II phase oscillators with analytical phase response curve (PRC) Z(ϕ) = −sin(ϕ): where subtracting these two relations gives the evolution equation for the relative phase of the two neurons: 12 21 where χ = ϕ 2 − ϕ 1 is the phase lag, and Ω = ω 2 − ω 1 . Using the trigonometric identity sin(x ± y) = sinx cosy ± cosx siny, Eq. (6) can be rearranged to Eq. (1).
Dynamical Analysis of the Joint Phase Model. In general, the fixed point of the phase lag χ ⁎ i of the Eq. (6) for type-I neuronal oscillator with analytical PRC function Z(ψ ± χ) = 1 − cos(ψ ± χ), is given by:  where χ ⁎ 1 is the inphase firing solution and χ ⁎ 2 belongs to antiphase state. Given the synaptic strengths, only one of these fixed points are stable in a given Ω and delay time ψ. The fixed points of type-II neuronal oscillator with analytical PRC function ψ χ ψ χ ± = − ± Z( ) sin( ) can be derived similarly: Eqs (7) and (8) show that the fixed points of both type-I and type-II neuronal oscillations are self-consistent in the presence of intrinsic frequency mismatch Ω between the two oscillators. In this case, the stable fixed point χ ⁎ numerically using any root-finding scheme. However, without loss of generality, one can assume that the intrinsic frequency mismatch between the two oscillators is negligible, i.e. Ω  0. Ignoring Ω, the type-I fixed point is still self-consistent, but the type-II fixed point can be simplified to:  where Γ is a positive expression that reflects the relative synaptic strength.
Network Model. A fully connected recurrent excitatory network with N = 200 type-II neuronal phase oscillators were simulated in the high-frequency regime with ν = 80 Hz. In the network simulation, the initial mean coupling is fixed at g (0) = 0.5. The elements of the adjacency matrix A are multiplied by the corresponding array in the synaptic strength matrix G in order to establish a weighted adjacency matrix A G . Note that the main diagonal arrays of the weighted adjacency matrix are zero since there is no self-loop in the network. Phase oscillators obey Eq. (1) and the excitatory synapses are modified by the pair-based STDP profile according to Eq. (2). The phases of the oscillators are initially uniformly distributed between 0 and π. The dendritic propagation delay is fixed at τ d = 0.5 ms. STDP parameters are A + = A − = 0.005, and τ + = τ − = 20 ms. We also define an order parameter r(t) for the network of phase oscillators ranging between 0 and 1 that measures the degree to which the system is synchronized 71 : Counting Loops. A bidirectional connection corresponds to a closed loop of length 2 in a network of neurons. In order to measure the number of such closed loops, a directed graph is constructed. Transformation of the synaptic strength matrix G, into a directed graph is performed by considering a threshold h = 0.2. Assuming that there are no self-loops, i.e. = g 0 ii , the corresponding value in the graph adjacency matrix M of the resultant directed graph is assigned to 1 when the synaptic strength is greater than h, and zero otherwise. The number of closed loops of length 2 in the graph adjacency matrix M is then 21,23 : where Tr denotes the matrix trace. In order to perform a better comparison, L 2 is normalized to the total number of possible loops of the same length i.e. N(N − 1)/2, ignoring self-loops, where N denotes the number of the neuronal phase oscillators or nodes in the network. Therefore, the resulted L 2 is a normalized number between 0 and 1.

Asymmetry Index.
To quantify the proportion of unidirectional connections and discriminate asymmetric modification of the synapses from symmetric changes (either symmetric potentiation or symmetric depression), we used another measure which has been introduced by Bayati and Valizadeh 17 . We first defined the synaptic cost of the network as the sum of all synaptic strengths in the weighted adjacency matrix A G , i.e. = ∑ S a g i j ij ij , , where a ij and g ij are the arrays of adjacency A and synaptic strength G matrices, respectively. The asymmetry level of the connections between two neurons can be measured by calculating the quantity C ij = −C ji = a ij g ij − a ji g ji . The asymmetry index of the network can then be defined as = ∑ | | > C S C (1/ ) i j ij net . By definition, the network asymmetry index C net is scaled in the range [0,1]. In the case of a fully symmetric network, C net = 0, whereas in a fully asymmetric network, C net = 1.