We utilize the network model from Ahn and Rubchinsky (2017) except that the synapses are plastic in this study. The model is described below.

The neurons are modeled using a two-dimensional conductance-based model of a Hodgkin-Huxley type that is mathematically equivalent to the Morris-Lecar model (Izhikevich, 2007; Ermentrout and Terman, 2010). The sodium conductance is assumed to activate instantaneously and to have no inactivation, while the potassium conductance is controlled by its gating variable and so varies dynamically.

Here v is the neuron's transmembrane potential and w is the gating variable for the potassium current. The synaptic current between neurons, I_syn, is given below and I_app is a constant input current to each neuron to control the frequency of spiking. The sodium, potassium, and leak currents are:

g_Na, g_K, and g_L are the maximal conductances for the sodium, potassium, and leak currents, respectively. The steady-state values for the gating variables of the sodium and potassium currents are:

The voltage-dependent activation time constant of the potassium current is:

All synapses are excitatory, and the synaptic current to neuron i is given by:

Where g_syn is the maximal conductance of the synapse (i.e., the synaptic strength), and s_j is the synaptic variable for neuron j and the summation is taken over all neurons that are connected to the i-th neuron. The synaptic variable s is governed by:

H_∞ is a sigmoidal function whose input is the presynaptic neuronal voltage:

The values of cellular and synaptic parameters are the same as used in Ahn and Rubchinsky (2017): g_Na = 1, g_K = 3.1, g_L = 0.5, v_Na = 1, v_K = −0.7, v_L = −0.4, v_m1 = −0.01, v_m2 = 0.15, v_w1 = 0.08, β = 0.145, I_app = 0.045, ε₁ = 0.02, ε₂ = 1.2ε₁, v_syn = 0.5, ∝_s = 5, β_s = 0.2, θ_v = 0.0, σ_s = 0.2.

STDP modeling follows (Zhigulin et al., 2003). If neuron i spikes at time t_i and neuron j spikes at time t_j, then the strength of the synapse from neuron i to neuron j is additively updated by the amount

where Δt = t_j − t_i. The synaptic conductance from neuron j to neuron i is simultaneously updated by an equal, but opposite, amount. While the additive update rule does not necessary need to be symmetric (as it is here), there is experimental evidence supporting the nature of the update, see for example (Zhang et al., 1998; Feldman, 2012). We varied the values of our plastic parameters, in particular A ∈ [0.0001, 0.01], k ∈ [0.01, 50]. The synaptic conductance is bounded below by zero.

The system of differential equations was solved numerically in Python using the built-in odeint function from the SciPy module (v.1.4.1). This function implements either the Adams method or a backward differentiation formula (BDF) method depending on the stiffness of the problem. The solution was reported at multiples of the time step dt = 0.1 (assuming the time units are milliseconds), however the function uses an adaptive step size and there was no lower bound on the length of the intermediate time steps that may be used (similarly, there was no upper bound restriction on the number of intermediate steps that were taken). The absolute and relative tolerances for the method were kept at the default value of 1.49 × 10⁻⁸. While the solution depends on the initial conditions, its statistical properties (such as the firing rate, synchrony pattern characteristics etc.) do not. The system was solved on the time interval [0, 25000], the first 20% of the time-series was removed from analysis. To implement plasticity, the integration was paused after each time step and, if necessary, the synaptic strength was updated. Specifically, the voltage threshold to define an action potential was set at 0.2.

The time-series analysis of synchronized dynamics in the network follows that of Ahn et al. (2011) and Ahn and Rubchinsky (2017) and is similar to the analysis of the temporal patterns of neural synchrony in the experimental studies mentioned in the Introduction. We will briefly describe this analysis here.

The phase, φ(t), of a neuron is defined as

where (ŵ,v^) is a point selected inside the neuron's limit cycle in the (w, v) – plane. The synchronization strength is computed as

where Δφ(t_j) = φ₁(t_j) − φ₂(t_j) is the difference of the phases of neurons 1 and 2 at time t_j. N is the number of data points. The value of γ ranges from 0 to 1, which represent a complete lack of synchrony and perfect phase synchrony, respectively.

If there is some degree of phase locking present, then there is a synchronized state, i.e., a preferred value of the phase difference Δφ. For each cycle of oscillation one can check if the actual phase difference is close to this preferred value or not. Note that the index γ only represents an average value of phase-locking over the interval [t₁, t_N], however to describe the patterning of synchrony one needs to look at the transitions to and from a synchronized state on much shorter timescales. This is done as follows.

When φ₁ increases past zero, say at time t_j,i, then φ₂(t_j,i) is recorded. This generates a sequence of numbers {φ2(tj,i)}i=1M. Due to the presence of some synchrony, there is a clustering about some phase value, say φ₀. This is taken as the preferred phase value, and if φ₂(t_j,i) = φ_i, for 1 ≤ i ≤ M, differs from it by more than π2 then the neurons are desynchronized, otherwise they are synchronized. The choice of π2 is not only convenient (it partitions the (φ_i, φ_i+1) space into quadrants) but was also used in the experimental studies described in the section Introduction.

The length of a desynchronization event is defined as the number of consecutive times the system spends in the desynchronized states. In other words, the length of desynchronization is the length of the time interval the system is away from the synchronized state (as defined above); this length is measured not in the absolute time units, but in the number of cycles of oscillations (in line with the experimental studies mentioned in the Introduction). The lengths of all desynchronization events are recorded and the distribution of durations is reconstructed. The mode of this distribution is used as a characteristic of the temporal patterning of synchronized dynamics. For later reference, a “mode n” system means that the mode of all lengths of desynchronization events for that particular system is n. Thus, a mode 1 system (n = 1 case) is the system with synchronized dynamics interrupted by predominantly short desynchronization intervals. The larger n is, the more prominent the tendency for long desynchronizations is. This does not necessarily affect the overall synchrony strength, because it depends not only on the duration of desynchronizations, but also on their number. The mode is used to characterize the durations because experimental studies used the mode for this purpose.

An illustration of different desynchronization durations and dynamics with different modes of desynchronizations is provided in Figure 1. Voltages and distributions of desynchronization durations for mode 1 dynamics are in the left column, the ones for mode 2 dynamics are in the right column. The synchronization is not perfect and synchronized dynamics (phase difference is close to the preferred one) are interspersed with desynchronized intervals. Note that the preferred phase difference is not necessarily zero so that the zero lag state is not necessarily a synchronized state.

Finally, we would like to reiterate that in this approach the time is measured in terms of cycles of oscillations of the neural activity, not in absolute time units. This allows one to compare the properties of variability of synchrony of brain rhythms with different frequencies.

The phase-locking strength index γ was observed to be usually about 0.2–0.3 in this study (even after STDP adjustments). These are moderate values, comparable with experimental results (in particular with the results reported in the studies references in the section Introduction). With this moderate synchrony strength, synchronization effects are hard to see by the naked eye, however, the quantitative time-series analysis techniques are able to quantify the synchronized dynamics and its properties including the temporal patterning of weakly synchronous dynamics.

Let us mention here that ε∝1τ and the maximum value of τ(v) is 1ε. Hence as ε is increased, the value of τ(v) is decreased across its entire domain as it is a unimodal function. This in turn accelerates the activation of potassium current because dwdt∝1τ(v). From Ahn and Rubchinsky (2017) we know that smaller values of ε promote shorter desynchronization events.

For ε = 0.05, the non-plastic system is mode 1. This means the synchronized dynamics has the following property. As the system is exhibiting partially synchronized dynamics, it will be either close in the synchronized state or away from synchronized state, the latter is termed desynchronization. The desynchronized interval length (measured in the number of cycles of oscillations) varies in time. We obtain the distribution of the desynchronization durations from numerical simulation and find the mode of this distribution. If this mode equals one cycle of oscillation, then the system is mode 1 (see section Methods for a more detailed explanation). Mode 1 means the desynchronizations are predominantly short.

Now the non-plastic system is changed to include STDP. The changes in the temporal patterning of synchronization dynamics are illustrated in Figure 4. Figure 4A is a diagram of the mode of the desynchronization durations in the space of plasticity parameters, A and k. The plasticity effects are negligible across the top (very large k implies a quick decay of the change in synaptic strength), and especially in the upper left corner (large k and a small amplitude A). In these areas the values of the plasticity parameters are such that the magnitude of the update, Δg_syn, is negligible (the average update is usually in the interval [0.0, 10⁻⁵], on the larger end this corresponds to about 0.2% of the initial value of g_syn). Hence, the plastic system continues to be mode 1 in these areas.

The rest of the parameter space, in particular the central region, displays a high proportion of mode 1 dynamics as well. In these areas plasticity is not negligible, as the synaptic strength can vary to a substantial degree. However, even in the presence of STDP, mode 1 dynamics persist. For the diagram in Figure 4A, about 85% of the parameter space points correspond to mode 1 systems.

To illustrate the effect of plasticity on a distribution of desynchronization durations, refer to Figures 4B–D. Plasticity effects increase from left to right. The distribution of durations changes: at a weak level of plasticity the durations are exclusively length one, while at a stronger level of plasticity some longer durations are observed. Yet the preponderance of length one desynchronization durations is preserved.

Now let us look at the effect of plasticity on the dynamics in systems with a mode larger than one. We consider ε = 0.15. The non-plastic system is mode 2 (the synchronization index γ is virtually unchanged from that of ε = 0.05, although the frequency of oscillations increases by several times, Ahn and Rubchinsky, 2017). Mode 2 means the desynchronizations tend to be longer than those of the mode 1 case.

Figure 5 shows the effect of STDP on the system that is mode 2 in the non-plastic case. As explained earlier, the plasticity effects are negligible across the top of Figure 5A, and especially in the upper left corner. We note that this region of the parameter space exhibits mode 2 dynamics (as expected). However, throughout the entire parameter space it is seen that a majority of parameter values correspond to mode 1 systems (the large central region in Figure 5A). Overall, about 20% of the parameter space points stay mode 2, while over 65% exhibit mode 1 dynamics (and less than 15% correspond to larger than mode 2 systems).

To illustrate the effect of plasticity on a distribution of desynchronization durations, refer to Figures 5B–D. Plasticity effects increase from left to right. Here we see that the introduction of weak plasticity can be sufficient to shift the system from mode 2 to mode 1 (Figure 5C). This means desynchronizations tend to become shorter in the plastic case. At stronger levels of plasticity (Figure 5D), the distribution widens, however the vast majority of desynchronization events remain length one.

Overall, we have seen that mode 1 dynamics are generally preserved when STDP is introduced to a non-plastic mode 1 system. When STDP is introduced to a non-plastic mode 2 system, the dynamics largely shifts from mode 2 to mode 1. The same was found with other non-plastic systems exhibiting higher modes: the introduction of STDP generally shifts the mode of the system down to one. Finally, we would like to note that there are several points in the parameter space (see Figures 4A, 5A) that have very large modes. For example, in Figure 4A when A = 0.0006, k = 0.01, the resulting system is mode 38 (i.e., most common desynchronizations are very long). Generally, these cases have a wide distribution of desynchronization durations. Therefore, while these systems have a large mode, the mode does not present a strong tendency in the distribution. Nevertheless, these situations are relatively rarely found.

The parameter β changes the widths of the voltage-dependent time-constant of activation τ(v) and the width of the steady-state activation function w_∞(v) for potassium current. In particular, as β is decreased, the slope at the half-height of w_∞(v) is increased, and this decreases the width of the step (w_∞(v) is a sigmoidal function). Similarly, for τ(v), a decrease in β decreases the width of the function around the peak. This causes an advancement in the activation of the potassium current.

A larger value of β promotes shorter desynchronization durations (Ahn and Rubchinsky, 2017). For β = 0.124, the non-plastic system is mode 1. The effect of STDP on this system is presented in Figure 6. Across the top and in the upper left corner of Figure 6A we see that virtually every point corresponds to a mode 1 system, as expected. Indeed, a substantial portion of the entire parameter space displays mode 1 dynamics; about 80% of the parameter space studied.

To illustrate the effect of plasticity on a distribution of desynchronization durations, refer to Figures 6B–D. Plasticity effects increase from left to right. The introduction of plasticity has a minimal effect on the distribution; there is very little change visibly. Indeed, the proportion of desynchronization durations of length one increases with plasticity.

Decreasing β increases the mode of a system. If β = 0.091, the non-plastic system is mode 2. With the introduction of very weak plasticity (across the top and the upper left corner of Figure 7A) we see that the dynamics are relatively unchanged, i.e., the mode of most systems remains two. However, if plasticity is not very weak, the dynamics shift to mode 1 in a significant portion of the parameter space. The effect is not as substantial as in the previous section, but about 35% of parameter space becomes mode 1 (about 45% remains mode 2, i.e., the mode is unchanged).

To illustrate the effect of plasticity on a distribution of desynchronization durations, refer to Figures 7B–D. Plasticity effects increase from left to right. We see that the vast majority of desynchronization durations become length one as plasticity becomes stronger.

The parameter v_w1 affects a horizontal translation in w_∞(v) and τ(v). Increasing v_w1 shifts both curves to the right, i.e., toward higher voltages; this results in a potassium current that activates faster.

Smaller values of v_w1 result in short desynchronization durations (Ahn and Rubchinsky, 2017). For v_w1 = 0.102, the non-plastic system is mode 1. The effect of STDP on this system is presented in Figure 8. We see that mode 1 dynamics is observed not only for the weak plasticity region (top and upper left corner of Figure 8A), but for most of the parameter space (about 85% of the parameter space studied).

To illustrate the effect of plasticity on a distribution of desynchronization durations, refer to Figures 8B–D. Plasticity effects increase from left to right. We see that as plasticity increases to a higher level, the prevalence of mode 1 is unchanged.

Varying v_w1 to larger values leads to shorter desynchronization durations becoming less prevalent. For v_w1 = 0.161, the non-plastic system is mode 2. The effect of STDP is presented in Figure 9. When plasticity is added we see that the dynamics are similar to the non-plastic case when plasticity is weak enough (top and upper left corner of Figure 9A). However, when the plasticity effects are moderate, the system exhibits mode 1 dynamics frequently (central region of Figure 9A). For the domain of parameter space studied, the majority of points (about 45%) correspond to mode 1 systems, the rest are either mode 2 (about 40%) or higher.

To illustrate the effect of plasticity on a distribution of desynchronization durations, refer to Figures 9B–D. Plasticity effects increase from left to right. We see that the mode of the system shifts down from two to one as plasticity becomes stronger.

Varying either ε, β, or v_w1 may affect not only the durations of the desynchronizations, but also synchronization strength and the frequency of activity in the system. To control desynchronization durations while keeping both spiking frequency and synchronization strength near constant in a non-plastic system, one can consider the parameter β and separate it into two independent parameters, β_τ and β_w. As a result, the lengths of desynchronization events are almost independent of the frequency and synchrony strength (Ahn and Rubchinsky, 2017).

Smaller β_w and larger β_τ result in shorter desynchronization durations (Ahn and Rubchinsky, 2017). For β_w = 0.098, β_τ = 0.079, the non-plastic system is mode 1. Figure 10 illustrates the impact of STDP on this system. Mode 1 dynamics is observed not only for the weak plasticity region (top and upper left corner of Figure 10A), but for the majority of the parameter space (about 60% of the parameter space studied).

To illustrate the effect of plasticity on a distribution of desynchronization durations, refer to Figures 10B–D. Plasticity effects increase from left to right. We see that as plasticity progresses to a moderate level, the proportion of short desynchronizations stays largely unchanged. In particular, the system is still mode 1.

If β_w = 0.115, β_τ = 0.071, the non-plastic system is mode 2. Figure 11 illustrates the impact of STDP on this system. With the addition of plasticity, we see that the system is largely mode 2 if the plasticity is weak (top and upper left corner of Figure 11A). However, stronger plasticity shifts the dynamics to mode 1 for a substantial portion of the parameter space (about 55% of points considered).

To illustrate the effect of plasticity on a distribution of desynchronization durations, refer to Figures 11B–D. Plasticity effects increase from left to right. We see that the distribution is largely unchanged for very weak plasticity, but as plasticity increases, the system becomes mode 1.