Synaptic reorganization of synchronized neuronal networks with synaptic weight and structural plasticity

Abnormally strong neural synchronization may impair brain function, as observed in several brain disorders. We computationally study how neuronal dynamics, synaptic weights, and network structure co-emerge, in particular, during (de)synchronization processes and how they are affected by external perturbation. To investigate the impact of different types of plasticity mechanisms, we combine a network of excitatory integrate-and-fire neurons with different synaptic weight and/or structural plasticity mechanisms: (i) only spike-timing-dependent plasticity (STDP), (ii) only homeostatic structural plasticity (hSP), i.e., without weight-dependent pruning and without STDP, (iii) a combination of STDP and hSP, i.e., without weight-dependent pruning, and (iv) a combination of STDP and structural plasticity (SP) that includes hSP and weight-dependent pruning. To accommodate the diverse time scales of neuronal firing, STDP, and SP, we introduce a simple stochastic SP model, enabling detailed numerical analyses. With tools from network theory, we reveal that structural reorganization may remarkably enhance the network’s level of synchrony. When weaker contacts are preferentially eliminated by weight-dependent pruning, synchrony is achieved with significantly sparser connections than in randomly structured networks in the STDP-only model. In particular, the strengthening of contacts from neurons with higher natural firing rates to those with lower rates and the weakening of contacts in the opposite direction, followed by selective removal of weak contacts, allows for strong synchrony with fewer connections. This activity-led network reorganization results in the emergence of degree-frequency, degree-degree correlations, and a mixture of degree assortativity. We compare the stimulation-induced desynchronization of synchronized states in the STDP-only model (i) with the desynchronization of models (iii) and (iv). The latter require stimuli of significantly higher intensity to achieve long-term desynchronization. These findings may inform future pre-clinical and clinical studies with invasive or non-invasive stimulus modalities aiming at inducing long-lasting relief of symptoms, e.g., in Parkinson’s disease.

In the original model [1], a neuron's average activity is depicted by its intracellular calcium concentration in accordance with experimental findings.However, the time-averaged firing rate could also represent the average neuronal activity since the calcium concentration and the firing rate are directly proportional [4].We use the time-averaged firing rate here, consistent with the stochastic SP model in the main text.
We simplified the BvOSP model as follows.The birth and death of synaptic elements of neuron i depend on its firing rate, f i , and the homeostatic set-point firing rate (target firing rate), f T .The synaptic elements are created when the activity of a neuron falls below f T and are deleted if it goes above f T .We assume a linear dependence of the change in the number of synaptic elements on f i .
where z represents both axonal and dendritic elements, i.e., z ∈ {Ax, Dn}.ν sp is the sprouting rate of the synaptic elements.We assume the same sprouting rate for both axonal and dendritic elements.At any given time, a synaptic element may be either bounded (as part of an existing synaptic contact) or vacant (available to form contacts).The detailed algorithm of addition and pruning of synaptic contacts is presented in Refs.[1][2][3].Briefly, dz i is calculated for each neuron until the following SP iteration and rounded to the smaller integer number.If dz i > 0, that many elements are created and added to the pool of vacant elements, while if dz i < 0, that many previously existing elements of neuron i are randomly selected and deleted, irrespective of their state (bounded or vacant).If the deleted element belonged to a synaptic contact, 1/12 the corresponding synaptic contact is pruned and the counterpart element becomes vacant and available for contact formation.The number of vacant elements decays spontaneously with a time constant, τ vac [1], The addition process is carried out by simply collecting the vacant synaptic elements and randomly assigning them to their counterparts, which causes a given postsynaptic (presynaptic) neuron to be more likely to connect to a presynaptic (postsynaptic) neuron with a larger number of vacant axonal (dendritic) elements [1].The contacts are formed if the distance-dependent probability permits.We use the same exponentially decaying distance dependence here as for the stochastic SP in the main text for the formation of a contact.Besides the conceptual difference, a key difference between BvOSP and the stochastic SP model is that BvOSP determines the number of outgoing contacts from presynaptic neurons based on their activity level via activity-dependent axonal bouton generation, while the stochastic SP model does not control the number of outgoing contacts.
BvOSP models only the homeostatic SP, meaning it does not include synaptic weight-dependent pruning.
We introduce the weight-dependent pruning of synaptic contacts, in which both axonal and dendritic elements engaged in that contact get deleted, to make it comparable to the stochastic SP model.A maximum of M contacts can be built between a pair of pre-and post-synaptic partners, i.e., k = 1, 2, .., M, and each contact between a pair can have a different weight at any given time.The change in the synaptic weights is governed by Eq. 7 in the main text.Only the contacts that become weak enough to get eliminated are pruned with the weight-dependent pruning probability (Eq.13 in the main text), while the others between the same pair remain intact.
Here we replace the homeostatic addition and pruning of our stochastic SP model introduced in the main text with BvOSP to compare the emerging network states and structure resulting from the two distinct models.As done in the main text using the stochastic SP model, we consider networks with 1. homeostatic SP (hSP) alone, i.e., without weight-dependent pruning (P w = 0) and without STDP, labeled as hSP-only; 2. a combination of STDP and hSP, i.e., without weight-dependent pruning (P w = 0), labeled as STDP+hSP ; 3. a combination of STDP and SP that includes hSP and weight-dependent pruning (P w ̸ = 0), labeled as

STDP+SP,
The BvOSP model is implemented as follows.The network is either initialized with random graph 2/12 connectivity or completely unconnected.dz i is calculated at every time-step and allowed to change over the window of duration ∆T sp during which no SP updates are executed and the network is allowed to settle in a steady state either with or without STDP, determined by Eq. 19 in the main text.The sprouting rate, ν sp = 10 −5 ms −1 and τ vac = 10∆T sp , comparable to previous studies [1][2][3][4].The addition of a synaptic contact from a neuron j to i increases the adjacency matrix element A i,j by 1 while the pruning of that reduces A i,j by 1 with the upper limit at M = 10 and lower at 0. The newly added contacts are given a random weight between 0 and 0.2 unless stated otherwise.

Networks with hSP-only
With BvOSP, we obtain similar results as with the stochastic SP method in the main text.First, we develop networks from completely unconnected neurons employing hSP-only with a given target firing rate, f T .We Neurons with a natural firing rate lower than f T develop more dendritic and axonal elements, leading to larger values of in-and out-NDDs, respectively.We observe an identical decrease of both in-and out-NDDs with an increase in the natural firing rates of the neurons in

Networks with combinations of structural and weight plasticity
We consider synchronized random networks with STDP-only and apply structural updates using STDP+SP (P w = 0.01), as outlined in the "Parameters and implementation" section in the main text, to investigate the effect of SP on such networks.We consider three random networks with initial average NDDs, β 0 = 0.075, 0.11,  using BvOSP that the peak at 0 disappears in the presence of weight-dependent pruning but not for P w = 0.
The number of strong contacts in the steady states for all three cases remains high.

Statistics of synaptic contacts between a pair of pre-and post-synaptic partners
We restrict the number of contacts from a presynaptic neuron to each of its postsynaptic partners with BvOSP to a maximum of 10, as previously suggested [5,6].established from a presynaptic neuron to each of its postsynaptic partners in three cases: networks developed from unconnected neurons using hSP-only for given target firing rates, f T , synchronized (initially) random networks with different initial NDDs, β 0 , that evolved with STDP+SP (P w = 0.01), and those that evolved with either P w = 0.01 or P w = 0.The counts remain much smaller than 10, and thus, may be unaffected by the imposed maximum, M, as previously observed [6].Previous studies that used their specially designed SP methods successfully reproduced experimentally observed bimodal distributions [5][6][7], where one peak is observed at 0 and the other tends to lie between 3 and 8 [5,8,9].However, we observed a consistently decreasing distribution with the increase in the number of contacts for all the cases we considered: networks generated from completely unconnected neurons using hSP-only as well as networks that were initially random but evolved with combinations of STDP and SP or hSP.It is because the multiple contact formation in our study corresponds to the independent synaptic contact formation via a random conversion of potential synaptic contacts to actual synaptic contacts, as previously discussed in Ref. [5].The distribution of potential synaptic contacts shows a consistent decrease, and so does the actual contact distribution in the present study.A bimodal distribution can be reproduced when the formation of new synaptic contacts between a given pair of neurons is assumed to depend on the number and state of the existing contacts, the specific interaction of neuron activity level, synaptic weight, and rewiring, and/or the presence of other neurons that may compete with the presynaptic partner [5][6][7]10].

All panels in
Here we considered no such dependencies and BvOSP inherently does not take into account the competition between the possible presynaptic partners [1].
Simulation results with networks of 32 × 32 neurons (N = 1024) In this section, we use the stochastic SP method instead of the BvOSP to show that the results presented in the main text with networks of 400 neurons hold for larger networks.We present an example where we develop networks using STDP+SP (P w = 0.01) and hSP-only to compare their transition to synchronized states.With STDP+SP, the new contacts added at each SP iteration (structural update) are given a random initial weight ∈ [0, 0.2] and are subjected to STDP allowing those to get potentiated or depressed.Thus, the 9/12 neurons can adjust their firing rates by adding and losing contacts besides changing the synaptic weights.
With hSP-only, the new contacts added at each SP iteration are given a random initial weight ∈ [0, 1].The target firing rate is set to 4.5 Hz.Since the target firing rate is higher than the natural firing rates of most of the neurons (mean 3 Hz and standard deviation 0.5 Hz), the initial structural updates involve only the addition of new contacts.The network may get synchronized when a sufficiently large number of contacts have developed with progressing SP iterations.This is indicated by an increase of the order parameter in consider f T = 3, 4, and 5 Hz.The newly added contacts are given a random weight between 0 and 1.The network-averaged firing rate increases with SP iterations as neurons develop synaptic contacts with other neurons to achieve f T , increasing the average NDD as shown in Fig A(a,c).Accordingly, a lower f T leads to a lower average NDD in the steady state, which may prevent the network from getting synchronized due to an insufficient number of contacts.Thus, a network with low f T settles in a desynchronized state, while those with larger values of f T may get partially or fully synchronized, as shown for f T = 4 and 5 Hz, respectively, in Fig A(b).A lower CV of firing rates accompanies a higher order parameter [Fig A(d)].
Fig B(a,b,c) for all values of f T considered since the number of axonal and dendritic elements follow the same dependence on f i and f T (Eq.1).The neurons with a natural firing rate above f T may remain unconnected, e.g., for f T = 3 and 4 Hz, as they may not develop any synaptic elements.Consequently, the joint probability distributions in Fig B(d,e,f) show a strong positive correlation between in-and out-NDDs.

Fig A.
Fig A. Time evolution of the state measures of networks developed with hSP-only starting from unconnected neurons.(a) shows the network-averaged firing rate, (b) the order parameter, (c) the network-averaged NDD, and (d) the CV of firing rates of the neurons.Parameter: sprouting rate, ν sp = 10 −5 ms −1 .

Fig C.
Fig C. Evolution of network-averaged measures with SP iterations.Colors correspond to the three initial values of average NDD (see legend).(a) shows the network-averaged firing rate, (b) the order parameter, (c) the network-averaged NDD, and (d) the CV of firing rates of the neurons.The parameters are: f T = 4.5 Hz, P w = 0.01, ⟨W ⟩(0) = 0.8.

Fig E.
Fig E. Dependence of average incoming and outgoing synaptic weight on the natural firing rate of neurons for a random network with STDP-only (a) and the distributions of the synaptic weight for random network with STDP-only, STDP+hSP (P w = 0), and STDP+SP (P w = 0.01) (b).Other parameters are the same as for Fig D.

Eq. 1
dictates that the lower the natural firing rate compared to the target rate, the higher the numbers of axonal and dendritic elements of the neuron.Axonal and dendritic elements are respectively responsible for outgoing and incoming synaptic contacts of a given neuron.With STDP+SP (P w = 0.01), out-NDD becomes low and in-NDD remains high for neurons with low natural firing rates because of the pruning of their weaker contacts in Fig F(a).The out-NDD increases and the in-NDD decreases with an increase in the 6/12 natural firing rate, following the synaptic weight dependence on the natural firing rate of the neurons shown in Fig E. With STDP+hSP (P w = 0), both in-and out-NDDs are high for neurons with low natural firing rates and show a minor drop with an increase in the natural firing rate of the neurons in Fig F(b) due to the absence of preferred pruning of weak contacts.Neurons with a natural firing rate close to the target rate have few to no synaptic elements.A small number of axonal elements leads to a low out-NDD, while that of dendritic elements results in a low in-NDD for both P w = 0.01 and P w = 0. Thus, we observe higher natural firing rate neurons with minimal in-and out-NDD in Fig F(a,b).

FigF
Fig F(c,d) show the joint probability distributions of the in-and out-NDDs of individual neurons for P w = 0.01 and P w = 0, respectively.For P w = 0.01, neurons with a high in-NDD have a low out-NDD, and most of the neurons with low in-NDD have a high out-NDD, resulting in a negative correlation between the two.For P w = 0, on the contrary, the in-and out-NDDs depend similarly on the natural firing rate of the neurons, and thus, assume similar values leading to a positive correlation between the two.

Fig F.
Fig F. Dependence of in-and out-NDD on the natural firing rates of the neurons for STDP+hSP (P w = 0), and STDP+SP (P w = 0.01) and their joint probability distribution, P (in-NDD,out-NDD) in the steady synchronized states.(a) and (b) show the dependence of the in-and out-NDD on the natural firing rates of the neurons for P w = 0.01 and P w = 0, respectively.(c) and (d) show the corresponding joint probability.Other parameters are the same as for Fig D.

Fig G
Fig G compares NDD distributions in steady synchronized states for P w = 0 and P w = 0.01.For P w = 0, both in-and out-NDD distributions have a similar shape and possess the same mode.The accumulation of Fig H shows the number of contacts (counts) Fig H(c) shows that the peaks at the non-zero numbers of contacts are higher for P w = 0.

Fig I(
Fig I(a) and a drop of the firing rate CV in Fig I(c).With hSP-only, the network may require a larger number of contacts to get synchronized because some of the strong contacts that are added may not support but rather oppose synchronization.However, for the network with STDP+SP, the synaptic weight change of the contacts may preclude the need to build more contacts, allowing for synchronization at a significantly earlier stage of network evolution and, thus, for smaller average NDD compared to the network with hSP-only.Combining panels (a) and (b) of Fig I and by excluding the SP iteration axis provides the dependencies of the order parameter on β in Fig I(d).Similarly, the dependencies of the firing rate CV on β in Fig I(e) are obtained by combining panels (c) and (b).We add the corresponding curves of the order parameter and the firing rate CV versus average NDD for random networks with STDP-only generated for every value of average NDD to compare the transition of networks developed from scratch to random networks with STDP-only.

Fig I(
Fig I(d) shows that the network with STDP+SP transitions to its synchronized state for significantly lower average NDD than networks with STDP-only and hSP-only.Fig I(b) shows that as the networks pass the onset of synchronization, the increase in average NDD remains smaller for the network with STDP+SP compared to that with hSP-only.This is further demonstrated in Fig I(d) and I(e), where the curves of order parameter and CV versus β cease to change at smaller values of β for networks with STDP+SP.Importantly, Fig I(e) shows that the plasticity of network structure results in synchronized states with much smaller values of the firing rate CV, indicating stronger frequency locking.

Fig I.
Fig I. Comparison of network evolution in the presence of different plasticity cases.The change in the order parameter, network-averaged NDD, and firing rate CV with SP iterations (a, b, and c, respectively).Order parameter and the CV of the firing rate versus the network-averaged NDD for the three plasticity cases (d and e).The asymmetry parameter, b, of STDP, is set to 1. Parameters: P h = 0.01, σ f = 0.5 Hz