Potential Mechanisms and Functions of Intermittent Neural Synchronization

Introduction

Synchrony of neural oscillations is believed to play important role in a variety of functions of the brain (e.g., Buzsáki and Draguhn, 2004; Colgin, 2011; Fell and Axmacher, 2011; Buzsáki and Schomburg, 2015; Fries, 2015; Harris and Gordon, 2015). Improperly organized (too excessive or too weak) synchrony is associated with several neurological and neuropsychiatric dysfunctions (e.g., Schnitzler and Gross, 2005; Uhlhaas and Singer, 2006, 2010; Oswal et al., 2013; Pittman-Polletta et al., 2015; Spellman and Gordon, 2015). However, the synchrony in cortical and subcortical circuits may not necessarily stay perfect for a prolong intervals of time (if it can be perfect at all). If synchrony is induced by transient stimulus, the transient character of this synchrony would be expected. But even in the idling dynamics of neural circuits of the brain prolonged perfect synchrony is rarely (if at all) reported.

This implies that for some intervals of time synchrony may be stronger, while for other intervals of time it may be weaker. The temporal patterns of synchrony may exhibit variations of synchrony strength yielding some average synchrony values. Few long intervals of desynchronized dynamics may be functionally different from many short desynchronized intervals, although the synchrony may be the same on the average.

Detection and quantification of the transient, varying, intermittent synchronization have been considered in the past (e.g., Hurtado et al., 2004; Le Van Quyen and Bragin, 2007). But these attempts were limited by the need in sufficiently long time-windows to obtain statistical significance, because synchronization is not an instantaneous phenomenon (Pikovsky et al., 2001). However, recent developments in time-series analysis (Ahn et al., 2011) allowed exploring the temporal patterning of synchrony on very short time-scales. The analysis of this fine temporal structure of synchronization is possible, because if some synchrony level is present on the average, then one can look at each cycle of oscillations and detect whether the signals are in the synchronous state or not.

These methods were used to study neural synchronization in several different systems: synchronization between single units and LFPs in the basal ganglia of Parkinsonian patients (Park et al., 2010; Ratnadurai-Giridharan et al., 2016), synchronization of EEG signals in healthy humans (Ahn and Rubchinsky, 2013), synchronization of LFPs between prefrontal and hippocampus circuits in normal rats and rats experiencing repetitive psychostimulant injections (Ahn et al., 2014). All these studies had one common feature: neural oscillations were observed to go out of synchrony frequently, but primarily for a short time interval. The observed synchrony level was reached by potentially very frequent, but short desynchronizations.

Since this short desynchronization dynamics was observed across different species, different conditions, and different signal types, it may be a universal feature of synchronized activity of neural systems. In this study we are providing a possible explanation for this apparent experimentally observed universality. We do so by looking for answers to two questions: what are the cellular or network mechanisms of this dynamics? What is its potential functional advantage?

We hypothesize that if this kind of dynamics is universal, it may be grounded in some general properties of neuronal excitability. In connection with this hypothesis, it is important to recall early insightful computational study (Somers and Kopell, 1993), which suggested that membrane conductances responsible for spiking help to speed up the establishment of synchrony. Here we will explore how experimentally observed short desynchronizations dynamics is defined by the kinetics of ionic channels, responsible for the generation of spikes. We also hypothesize that short desynchronization dynamics permits creation of synchronous states with weaker inputs. This may make neural systems more adaptable as they can easily create synchronous assemblies in response to synaptic or sensory inputs.

Since short desynchronization dynamics may be a generic phenomenon based on the properties of membrane channels, which are hard to alter in experiment, we use computational neuroscience techniques to study very simple conductance-based neuronal models. We alter the properties of conductances to explore their critical features for short desynchronization dynamics and investigate how coupled neurons may be efficiently entrained by external input. Models are subjected to the same kind of time-series analysis techniques as were used in earlier experimental studies. As a result, we reveal potential cellular basis of short desynchronization dynamics in the model and present its potential functional advantages.

Methods

Neuronal Model

We use a conductance-based modified Morris-Lecar neuronal model (Izhikevich, 2007; Ermentrout and Terman, 2010). We choose it because it is a simple (perhaps, the simplest) model that directly retains membrane conductances. Even though the original Morris-Lecar model includes calcium and potassium currents, it is equivalent to a reduced classical Hodgkin-Huxley sodium-potassium model (Izhikevich, 2007; Ermentrout and Terman, 2010). So, by studying this neural model, we model neurons with sodium-potassium spiking mechanism with fast sodium and delayed rectifier potassium currents.

We consider the model in the form:

$\begin{array}{llll}v\prime & =& -I_{Na}-I_K-I_L-I_{syn}+I_{app},& (1)\end{array}$

$\begin{array}{ll}w\prime =\left[w_{\infty }(v)-w\right]/\tau (v).& (2)\end{array}$

v is a transmembrane voltage and w is the gating variable of potassium current. I_Na = g_Na m_∞(v)(v − v_Na), I_K = g_Kw (v − v_K) and I_L = g_L(v − v_L) are sodium, potassium, and leak currents; I_app is a constant parameter and I_syn is a synaptic current (see below). g_Na, g_K, g_L are the maximal conductances for the Na⁺, K⁺, and the leak currents. The functions

$\begin{array}{ll}{m}_{\infty }\left(v\right)=\frac{1}{1+exp\left(-2\left(\frac{v-v_{m_1}}{v_{m_2}}\right)\right)},& (3)\end{array}$

$\begin{array}{ll}{w}_{\infty }\left(v\right)=\frac{1}{1+\text{exp}\left(-2\left(\frac{v-v_{w_1}}{\beta }\right)\right)},& (4)\end{array}$

$\begin{array}{ll}\tau \left(v\right)=\frac{1}{\epsilon }\frac{2}{exp\left(\frac{v-v_{w_1}}{2\beta }\right)+\text{exp}\left(\frac{-\left(v-v_{w_1}\right)}{2\beta }\right)}& (5)\end{array}$

are the steady-state activation functions of the gating variables of the Na⁺ and K⁺ currents, and the activation time-constant for K⁺ current. The functions m_∞(v) and w_∞(v) have sigmoid shapes while τ(v) has a unimodal shape. The term I_syn represents the synaptic current between cells.

We consider neurons connected with excitatory synapses adapted from Izhikevich (2007) and Ermentrout and Terman (2010). For a cell i, the synaptic current $I_{syn,i}=g_{syn}\left(v_i-v_{syn}\right){\displaystyle \sum _{j\ne i}}{s}_j$, where the sum is over those cells that send synaptic inputs to a cell i. The synaptic variable s is modeled by the first-order kinetic equation in the form:

$\begin{array}{ll}s\prime ={\alpha }_s(1-s)H_{\infty }(v-{\theta }_v)-{\beta }_{s}s,& (6)\end{array}$

where $H_{\infty }\left(x\right)=1/[1+exp\left(-\frac{x}{{\sigma }_s}\right)]$ is a sigmoidal function and v is the presynaptic neuron voltage (Izhikevich, 2007; Ermentrout and Terman, 2010). The parameter values are g_Na = 1, v_Na = 1, g_K = 3.1, v_K = −0.7, g_L = 0.5, v_L = −0.4, I_app = 0.045, v_m1 = −0.01, v_m2 = 0.15, β = 0.145, v_w1 = 0.08, ε = 0.02, g_syn = 0.005, v_syn = 0.5, α_s = 2, β_s = 0.2, θ_v = 0, σ_s = 0.2. We will further vary the values of ε, β, and v_w1 as will be described in the Results.

Synchronization Analysis

Phase analysis is frequently used to analyze synchronous neural activity of both continuous (LFP, EEG) and spiking signals (see, e.g., Lachaux et al., 1999; Hurtado et al., 2004; Le Van Quyen and Bragin, 2007). This analysis was used in experimental studies revealing prevalence of short desynchronization dynamics (Park et al., 2010; Ahn and Rubchinsky, 2013; Ahn et al., 2014; Ratnadurai-Giridharan et al., 2016). So we assume a very similar approach here. For a spiking activity, a phase of a neuron i is reconstructed by computing

$\begin{array}{ll}{\phi }_i\left(\text{t}\right)=\text{arctan}\left(\frac{v_i\left(t\right)-{\widehat{v}}_i}{w_i\left(t\right)-{ŵ}_i}\right),& (7)\end{array}$

where $({\widehat{v}}_i,{ŵ}_i)$ is a rest state of a neuron, set as a center of rotation in the (v_i, w_i)-plane. Then we consider an average synchronization index to measure the strength of the phase locking between two signals (Pikovsky et al., 2001; Hurtado et al., 2004):

$\begin{array}{ll}\gamma =\left|\right|\frac{1}{N}{\sum }_{j=1}^N{e}^{i\Delta \phi (t_j)}|{|}^2,& (8)\end{array}$

where Δφ(t_j) = φ₁(t_j) − φ₂(t_j) is the phase difference, the t_j are the sampling points, N is the number of data points to be considered, and ||.|| is the absolute value of a complex number.

This phase synchronization index γ varies from 0 (lack of synchrony) to 1 (perfect synchrony). It provides an average value of phase-locking. There may be cycles of oscillations, where phase difference is close to the average value of the phase difference (phase-locked, synchronized state) and where it is not close to it (desynchronized state).

To study the fine temporal structure of the dynamics of synchronization we use the methods recently developed in Park et al. (2010) and Ahn et al. (2011). Whenever φ₁ crossed the zero from negative to positive values, we recorded the value of φ₂, generating a set of consecutive phase values {ϕ_i}, i = 1, …, N. If the value of ϕ_i differs from the average value of ϕ_i by less than π/2 then the neurons are considered to be in a synchronized state, otherwise they are in the desynchronized state. We chose the value of the threshold to be π/2 because the experimental studies we discussed above used this value. The duration of desynchronization events is defined as the number of cycles of oscillations that the system spends in the desynchronized state minus one. The mode of the distributions of desynchronization durations is defined as the number with the highest probability of desynchronization durations.

We characterize the fine temporal structure of intermittent synchronization by quantifying the properties of distribution of desynchronization durations. We compute the relative frequencies (probabilities) of the durations of desynchronization events. This is similar to how the experimental data were characterized in the studies of the temporal patterns of synchrony (Park et al., 2010; Ahn and Rubchinsky, 2013; Ahn et al., 2014; Ratnadurai-Giridharan et al., 2016). We use the mode of the distribution of desynchronization durations and the probability to observe this mode, p_mode. If the mode of the desynchronization duration is short, but other desynchronizations (especially longer ones) are almost as frequent, then the dynamics is not necessarily dominated by short desynchronizations. However, if p_mode is close to one, then all other desynchronization durations are rare.

In our approach the duration of synchronization and desynchronization intervals is measured not in the absolute time units, but in cycles of oscillations, as was done in experimental studies. It makes easier to compare synchronization patterns between rhythms of different frequency. However, as we study the differences between different desynchronization durations in the modeling, we also compare the dynamics with the same frequencies of rhythms (see Results).

Results

We will study the dynamics of coupled model neurons as we vary parameters of potassium current. We do so by varying three different parameters: ε, β and v_w1 (see Equations 4 and 5), they all affect the effective value of activation time-constant τ(v) of potassium current. Larger values of τ delay activation of potassium current and promote characteristic shape of spike with very sharp rise of voltage, faster decay of voltage, and prolong interval between spikes. Lowering effective values of τ in the model (for example, by using larger values of ε) will lead a model neuron to generate less spiky and more quasi-sinusoidal profile of activity (Figure 1). By changing the values of ε, β, v_w1 we can study the model neurons exhibiting spiking activity like at Figure 1A as well as more sinusoidal activity like at Figure 1B (which is not necessarily very realistic, but will help in understanding the mechanisms and functions of physiological activity).

FIGURE 1

Figure 1. Numerical simulation of a voltage of an isolated neuron (Equations 1–5). Examples of spiking activity (ε = 0.001) (A) and quasi-sinusoidal activity (ε = 0.5) (B).

Kinetics of Voltage-Gated Potassium Channel and the Temporal Patterning of Synchronization

We consider a minimal neuronal network to exhibit synchronized dynamics: two neurons mutually connected with excitatory synapses (Figure 2). These two neurons satisfy (Equations 1–5) and the synapses are described by Equation (6). We consider two weakly coupled neurons with a small difference in firing rate (that is, the frequency of oscillations of voltage) due to slightly different ε. Note that since ε₁ ≠ ε₂ and the coupling strength g_syn = 0.005 is weak, two cells are not fully synchronized. Thus, the synchrony is intermittent rather than perfect.

FIGURE 2

Figure 2. Diagram of a minimal network of excitatory coupled neurons. We use ε₁ ≠ ε₂ (i.e., neurons have different firing rates) and the coupling strength g_syn is not very strong.

The Effect of the Peak Value of Activation Time-Constant

The magnitude of the voltage-dependent activation time-constant τ(v) of potassium current is inversely proportional to the parameter ε in such a way that the maximal value of τ(v) is 1/ε. We consider here how ε affects the durations of desynchronization events and accompanying changes in the average synchrony level and the mean frequency of spiking.

As the values of ε₁ and ε₂ increase, the fine temporal structure of synchronization changes as evident by the changes of the mode of the distribution of desynchronization durations (Figure 3A). Smaller values of ε promote short desynchronization intervals lasting for only one cycle of oscillations. On the contrary, the increase in ε leads to the increase of the mode of the distribution of desynchronization durations. That is, as ε increases, the most frequent desynchronization intervals are getting longer. Since the activation time-constant τ(v) is inversely proportional to ε, larger value of τ(v) (which promotes spike-like waveform in the model) promotes short desynchronization dynamics.

FIGURE 3

Figure 3. The effect of ε₁ when ε₂ = 1.2ε₁. (A) Mode value of the durations of desynchronization events (black line with black dots) and the corresponding probability to observe the mode value (gray line). The high (close to one) value of probability at mode indicates that the desynchronizations of corresponding duration are strongly prevalent. (B) Synchronization strength index γ. (C) The mean frequency of activities of both neurons. Since ε₂ = 1.2ε₁, neuron 2 has slightly higher frequency than the mean frequency while the neuron 1 has slightly lower frequency than the mean frequency.

The synchrony strength γ experiences only very small variations (Figure 3B). This indicates that the durations of desynchronizations may be independent of the synchrony strength. The same level of synchrony may be reached with numerous short desynchronizations or few long desynchronizations.

The mean frequency (firing rate) grows substantially (Figure 3C). This is expected because the growth rate of w is proportional to ε (Equation 2). As a result, while desynchronization intervals measured in cycles of oscillations are longer for larger ε, their durations in the absolute time are not necessarily growing. We will address this issue below.

Note that the probability of the dominant duration of desynchronization events p_mode (thin gray line without dots in Figure 3A) is mostly close to 1 and always higher than 0.5. Thus, the desynchronization durations of the corresponding number of cycles is really dominant (because the sum of all probabilities of the durations of desynchronization events is 1).

The Effect of the Widths of Voltage-Dependence of the Activation Time-Constant τ(v)

In the model, the parameter β is related to the width of the steady state activation function w_∞(v) and range of activation constant τ(v), where τ(v) is substantially different from 0. The results presented below will show that it is mostly the widths of voltage-dependence of the activation time-constant τ(v) that matters for the properties of desynchronization durations. As β increases, the width of τ(v) increases. That is, the range of voltages, where the activation time-constant is different from 0, is getting larger. This may effectively bring τ closer to the maximal possible value (which is 1/ε in the model) for a larger range of voltage and thus for a longer time. Thus, similar to the decrease of ε, larger β will promote more “spiky” and less quasi-sinusoidal waveforms.

Figure 4 shows how the parameter β affects the synchronized dynamics of coupled neurons. Larger values of β promote shorter desynchronization episodes (Figure 4A). There is also an effect on the synchrony strength (Figure 4B) and the frequency (Figure 4C). Shorter desynchronizations correspond to the higher synchrony level. As β is changing, the frequency is changing. This may mitigate the short desynchronization phenomena if desynchronization duration is measured in absolute units of time instead of cycles of oscillations. Nevertheless, similar to the case considered above, a change in parameter that leads to larger values of activation time-constant τ (increase in β) promotes desynchronizations of shorter durations [as signified by the high value of the probability at the mode of distribution of desynchronization durations (gray thin line in Figure 4A)].

FIGURE 4

Figure 4. The effect of β when ε₁ = ε and ε₂ = 1.2ε. (A) Mode value of the durations of desynchronization events (black line with black dots) and the corresponding probability to observe the mode value (gray line). (B) Synchronization strength index γ. (C) The mean frequency of activities of both neurons.

The Effect of the Voltage of Half-Activation and Maximal Activation Time-Constant

We now consider the effect of v_w1 which is the midpoint of the steady state activation function w_∞(v), that is the voltage at which half of the channels open. The same parameter defines the voltage at which the activation time-constant τ(v) peaks. Increase in v_w1 shifts both curves w_∞(v) and τ(v) in the direction of higher voltages. The conductance will start to increase at higher voltage and at a later time, and will start to decrease at earlier time. Thus, larger value of v_w1 may be expected to have an effect analogous to decrease in τ.

The results of numerical simulation presented at the Figure 5 support this. Lower values of v_w1 promote shorter desynchronizations (Figure 5A). In this case, the synchronization strength is larger for short desynchronization dynamics (Figure 5B), but the frequency is almost constant (Figure 5C). Thus, the desynchronizations are short here not only if measured in the number of cycles (spikes), but also measured in absolute time units.

FIGURE 5

Figure 5. The effect of v_w1 when ε₁ = ε and ε₂ = 1.2ε. (A) Mode value of the durations of desynchronization events (black line with black dots) and the corresponding probability to observe the mode value (gray line). (B) Synchronization strength index γ. (C) The mean frequency of activities of both neurons.

Changing Desynchronization Durations Independently of Both Frequency and Synchrony Strength

The changes in desynchronization durations in the numerical experiments above are accompanied by changing either average synchrony strength or firing rate (or even both). Here we consider whether the desynchronization durations can vary independently of both synchrony strength and firing rate. To study this, we modify the Equations (4) and (5) so that the values of parameter β in equations for w_∞(v) and τ(v) are not identical. This means the half-activation voltage is different from the voltage at which τ(v) has maximal value:

$\begin{array}{ll}{w}_{\infty }\left(v\right)=\frac{1}{1+\text{exp}\left(-2\left(\frac{v-v_{w_1}}{{\beta }_w}\right)\right)},& (9)\end{array}$

$\begin{array}{ll}\tau \left(v\right)=\frac{1}{\epsilon }\frac{2}{exp\left(\frac{v-v_{w_1}}{2{\beta }_{\tau }}\right)+\text{exp}\left(\frac{-\left(v-v_{w_1}\right)}{2{\beta }_{\tau }}\right)}.& (10)\end{array}$

Smaller β_w makes the slope of the steady-state activation function w_∞(v) larger, while smaller β_τ makes the width of the constant function τ(v) smaller. We let β_w and β_τ be changing in opposite directions. As β_w decreases from 0.134, β_τ increases from 0.061 at a different rate (β_w = 0.134 − 0.001k and β_τ = 0.061+0.0005k, where k = 0, 1, …, 40.). For other parameters, we use ε₂ = 1.3ε₁, ε₁ = 0.03, I_app = 0.04, g_syn = 0.005, v_w1 = 0.07. These changes of β_w and β_τ are not necessarily biologically realistic, but they allow us to explore whether the changes of desynchronization durations must covary with the changes of average synchrony or firing rate.

Figure 6 shows that in this case the synchrony strength and the firing rate are almost constant while the mode of desynchronization durations changes drastically. In other words, simultaneous variations of the width of w_∞(v) and τ(v) vary the distribution of desynchronizations independently from synchrony strength and firing rate. Thus, the same level of synchrony strength may be supported either with many short desynchronizations or few long desynchronizations regardless of whether the durations of desynchronizations are measured in cycles of oscillations or in absolute time units.

FIGURE 6

Figure 6. Changing synchronization durations independently from synchrony strength and firing rate. (A) Mode value of the durations of desynchronization events (black line with black dots) and the corresponding probability to observe the mode value (gray line). (B) Synchronization strength index γ. (C) The mean frequency of activities of both neurons.

Short Desynchronization Dynamics and Synchronization Threshold

To study potential functional advantages of short desynchronizations dynamics, we will consider two mutually excitatory connected neurons (as before) receiving common synaptic input from a third neuron: neuron 3 excites neurons 1 and 2 through the excitatory synapses (but does not get any feedback, Figure 7).

FIGURE 7

Figure 7. Diagram of a minimal network of excitatory coupled neurons receiving common synaptic input. Neurons 1 and 2 have different firing rates. They are mutually coupled through excitatory synapses and receive synaptic input from neuron 3.

We consider two different versions of three-neuron networks in Figure 7. In the first version, the parameters are selected in such a way that when g_syn1 = 0 neurons 1 and 2 exhibit dynamics with mostly short desynchronizations. The second version exhibits partially synchronized dynamics with the most common desynchronization intervals lasting for 4 cycles of oscillations when g_syn1 = 0. In other words, we consider how two coupled neurons exhibiting either short desynchronization dynamics or longer desynchronization dynamics respond to the common synaptic input.

One network has β_w = 0.094 and β_τ = 0.081, this is the left end of the horizontal axis in Figure 6A. The mode of desynchronization durations is just 1 and we will call this network “cycle 1” network (short desynchronizations network). The other network has β_w = 0.134 and β_τ = 0.061, this is the right end of the horizontal axis in Figure 6A. The mode of desynchronization durations is 4 and we will call this network “cycle 4” network (longer desynchronizations network). It is important to note that both networks have almost the same synchrony strength (Figure 6B) and firing rate (Figure 6C). So, except the difference in desynchronization durations, the dynamics of two networks are similar. That is, they have the same synchrony level and the same period of oscillations in the absence of synaptic input from the neuron 3.

In the numerical experiments, ε₁ = 0.03 and ε₂ = 1.3ε₁, the same values as used in Figure 6. We consider four different values of the firing rate in the neuron 3: ${\epsilon }_3=\left\{0.5{\epsilon }_1,\frac{{\epsilon }_1+{\epsilon }_2}{2},1.5{\epsilon }_1,2{\epsilon }_1\right\}$. So, the firing rate in neuron 3 is either substantially lower than in neurons 1 and 2, equals to the average of firing rates of neurons 1 and 2, or is higher than firing rates in neurons 1 and 2. All other parameters of neurons 1, 2, and 3 are the same and fixed as those used in Figure 6.

Now let us consider these two networks as the common input to neurons 1 and 2 is getting stronger due to increase of g_syn1 from zero (while g_syn = 0.005 is fixed, that is, the coupling between neuron 1 and neuron 2 is relatively weak). As synaptic input from neuron 3 to neurons 1 and 2 is getting stronger, neurons 1 and 2 are becoming more synchronous and will eventually be in full synchrony with each other due to common synaptic input and mutual synaptic coupling.

We compute the synchrony index γ for “cycle 1” and “cycle 4” networks (γ^(1C) and γ^(4C) respectively) when increasing values of g_syn1. To study how differently these networks are synchronized, we consider the absolute and relative difference of synchronization indices γ^(1C) and γ^(4C) for different values of g_syn1. Figure 8 presents the averages of γ^(1C) − γ^(4C) (thick solid line) and $\frac{{\gamma }^{(1C)}-{\gamma }^{(4C)}}{{\gamma }^{(4C)}}$ (the insert figure). Both quantities indicate how much synchronization in “cycle 1” (short desynchronizations) network is stronger than synchronization in “cycle 4” network when they receive the same synaptic input from the neuron 3.

FIGURE 8

Figure 8. The synchrony difference γ^(1C) − γ^(4C) is plotted for different strength of synaptic input g_syn1, normalized synchrony difference ($\frac{{\gamma }^{(1C)}-{\gamma }^{(4C)}}{{\gamma }^{(4C)}}$) is presented at the inserts (γ^(1C) and γ^(4C) represent the synchrony index γ for “cycle 1” (short desynchronizations) and “cycle 4” (long desynchronizations) networks, respectively). Subplots (A–D) are for different values of the firing rate of incoming signal, corresponding to ${\epsilon }_3=\left\{0.5{\epsilon }_1,\frac{{\epsilon }_1+{\epsilon }_2}{2},1.5{\epsilon }_1,2{\epsilon }_1\right\}$ respectively.

When this input is weak (g_syn1 is small), γ^(1C) and γ^(4C) are close to each other. When g_syn1 is large, γ^(1C) and γ^(4C) are again close to each other because both networks are necessarily strongly synchronous due to strong input. However, for the values of g_syn1 between zero and synchronization threshold value, γ^(1C) − γ^(4C) is large and positive. So the networks exhibiting short desynchronizations dynamics in the absence of input (“cycle 1” networks) reach either the same synchrony levels or higher synchrony levels than long desynchronization (“cycle 4”) networks for the same strength of synaptic input g_syn1. This phenomenon is observed regardless of the firing rate in presynaptic neuron 3 (i.e., regardless of ε₃). Sometimes this difference in the synchronization strength is moderate, but sometimes it is quite substantial (see Figure 8). The synchrony index is bounded by one from above, so the magnitude of the phenomenon is more emphasized by observing the relative value of synchronization index difference (inserts in Figure 8).

We also measure the threshold value g_syn1 for two neruons to reach synchornized dynamics without desynchronization events (Figure 9). This does not imply complete synchrony, but implies only small deviations between the phases of two signals, so that it is small enough to have no desynchonization events. As can be seen in Figure 9, the computed synchrony thresholds for short desynchronization (“cycle 1”) network were lower than the synchrony thresholds for long desynchronization (“cycle 4”) network for all considered firing rates (all possible ε₃).

FIGURE 9

Figure 9. Threshold value of synaptic strength g_syn1 to reach synchornized dynamics without desynchronization events for different values of ε₃. Black squares represent the critical value of g_syn1 for short desynchronization (“cycle 1”) network and the gray circles represent the critical value of g_syn1 for long desynchronization (“cycle 4”) network.

The results presented in Figures 8, 9 indicate that with average synchrony level and mean firing rate being equal, neural systems with short desynchronization dynamics reach higher synchrony for the same synaptic input strength and need weaker inputs to be synchronized than neural systems with long desynchronization events.

Discussion

Cellular Mechanisms of Short Desynchronization Dynamics

Imperfect synchrony is widely observed in the activity of neural networks of the brain. New time-series analysis techniques showed that intervals of synchronous dynamics are interspersed between desynchronized episodes, and most desynchronized episodes are very short (see references in Introduction). This stereotyped fine temporal structure of neural synchronization is not an artifact of the analysis method because other types of patterning of synchronization are possible in non-neural coupled oscillators (Ahn et al., 2011; Rubchinsky et al., 2014).

The present study provides potential mechanisms for this type of temporal patterning of neural synchrony. We varied several parameters of potassium conductance and identified conditions leading to the intermittent neural synchrony with predominantly short desynchronization episodes similar to experimental ones. All these conditions (large peak value of activation time-constant, large width of dependence of activation time-constant on voltage, lower values of voltage for peak activation time-constant) lead to the relatively large values of the activation time-constant τ(v) in the right range of voltages. The large value of τ(v) leads to the delay in activation of potassium current, so that a sharp spike can be generated. And, as our results show, it promotes the short desynchronizations dynamics.

The results of the computational modeling also indicate that the distribution of desynchronization durations may be independent of the synchronization strength. The same synchrony strength may be achieved with desynchronizations of different durations. Moreover, our results regarding comparison of synchronization in networks exhibiting short desynchronizations and long desynchronizations are obtained for the case when not only average synchrony level is the same, but the period of oscillations (the firing rate) is the same. By appropriate adjustment of model parameters we dissociated the effects of frequency of oscillations and of average synchrony strength from the effects of fine temporal patterning of synchronized dynamics.

These model-based observations fit with experimental observations of the changes in the distribution of desynchronization durations in prefrontal cortex-hippocampal synchrony in behavioral sensitizations experiments (Ahn et al., 2014). In these experiments, the desynchronization durations were predominantly short and their distributions were altered after psychostimulant administration, while the average synchrony levels stayed the same.

The “spikiness” of oscillations of transmembrane voltage is a very generic property of many neurons, which relies on the fast activation of current with high reversal potential and slow activation of current with low reversal potential. Our results show that the same conditions that promote short desynchronization dynamics promote the characteristically sharp shape of an action potential. This may explain why very different neural systems exhibit short desynchronization dynamics, as we described in Introduction.

Limitations of the Modeling Approach

We use a very simple model of a neuron and very simple model of a network. There are many factors, which may affect synchronous dynamics of neural activity, yet they are not represented in the model. Other important factors, which affect neural synchrony, are different membrane currents and their properties (we have a model with just two conductances and consider only several parameters of one conductance) and the size of the network (we have a very small network). Heterogeneity of the networks is also important (we have a very minimal representation of heterogeneity). Synaptic plasticity is not incorporated in our model (and is the subject of the future research). Finally, noise may affect temporal patterns of synchrony, which is not considered in this study either.

However, even though these factors are not incorporated in the model (which captures only some very basic mechanisms of neural activity), the model is able to generate realistic synchrony patterns. So, the right way to interpret the modeling results is to see what these basic mechanisms are capable of. These modeling results suggest that these very basic neural mechanisms are capable of explaining the properties of experimentally observed intermittency of neural synchrony. As we discussed above, short desynchronization dynamics has been observed in several different neural systems. An ability of a minimal neural network considered here to describe the properties of the intermittent synchrony (which is common to all those systems) is probably an indicator that the general neural mechanisms built in the model are adequate to the considered phenomena.

Inhibition is playing an important role in neural synchronization, but is not considered in our model. The experimental data discussed here were collected from cortical and subcortical networks with excitatory and inhibitory synapses. It will be interesting to see how the intermittent patterns of synchrony are affected by inhibitory synapses.

We would also like to note that our earlier study with more advanced neural and network model (which included excitatory and inhibitory synapses) did provided a quantitatively adequate description of the short desynchronization dynamics at the beta-band oscillations in the basal ganglia in Parkinson's disease (Park et al., 2011). The present modeling study is not designed to provide a quantitative description of a specific experiment, but rather it provides a qualitative description of common aspects of neural synchrony in different neural systems.

Potential Functional Significance of Short Desynchronization Dynamics

Our computational results suggest one way of how short desynchronization dynamics can be beneficial for neural systems. With two important properties of dynamic (average synchrony strength and firing rate) being equal, neural systems with short desynchronizations are easier to synchronize with common synaptic input. We showed that the same strength of common synaptic input leads to larger synchrony level in short desynchronization system. In other words, short desynchronization dynamics allows reaching a pre-set synchrony level with weaker input. So, if a strong synchrony is needed, systems with short desynchronizations will reach the pre-set synchrony strength with weaker inputs compared to longer desynchronizations.

Given the functional importance of synchronization in many neural systems (see references in Introduction), short desynchronization dynamics may allow for efficient regulation of synchrony levels. While the same level of synchrony may potentially be achieved with few long desynchronization episodes as well as with many short desynchronization episodes, only short desynchronization dynamics is experimentally observed in the neural synchrony in the brains. Our modeling results suggest that this short desynchronizations dynamics is easier to control with synaptic input. Thus, very basic properties of delayed rectifier potassium current (its delayed activation) is likely to promote efficient formation and break-up of synchronized assemblies.

Author Contributions

LR conceived research; SA and LR designed research; SA performed numerical simulations, SA and LR analyzed and interpreted the results; SA and LR wrote the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Ahn, S., Park, C., and Rubchinsky, L. L. (2011). Detecting the temporal structure of intermittent phase locking. Phys. Rev. E 84:016201. doi: 10.1103/PhysRevE.84.016201

PubMed Abstract | CrossRef Full Text | Google Scholar

Ahn, S., and Rubchinsky, L. L. (2013). Short desynchronization episodes prevail in synchronous dynamics of human brain rhythms. Chaos 23:013138. doi: 10.1063/1.4794793

PubMed Abstract | CrossRef Full Text | Google Scholar

Ahn, S., Rubchinsky, L. L., and Lapish, C. C. (2014). Dynamical reorganization of synchronous activity patterns in prefrontal cortex - hippocampus networks during behavioral sensitization. Cereb. Cortex 24, 2553–2561. doi: 10.1093/cercor/bht110

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsáki, G., and Draguhn, A. (2004). Neuronal oscillations in cortical networks. Science 304, 1926–1929. doi: 10.1126/science.1099745

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsáki, G., and Schomburg, E. W. (2015). What does gamma coherence tell us about inter-regional neural communication? Nat. Neurosci. 18, 484–489. doi: 10.1038/nn.3952

PubMed Abstract | CrossRef Full Text | Google Scholar

Colgin, L. L. (2011). Oscillations and hippocampal-prefrontal synchrony. Curr. Opin. Neurobiol. 21, 467–474. doi: 10.1016/j.conb.2011.04.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Ermentrout, G. B., and Terman, D. H. (2010). Mathematical Foundations of Neuroscience. New York, NY: Springer.

Google Scholar

Fell, J., and Axmacher, N. (2011). The role of phase synchronization in memory processes. Nat. Rev. Neurosci. 12, 105–118. doi: 10.1038/nrn2979

PubMed Abstract | CrossRef Full Text | Google Scholar

Fries, P. (2015). Rhythms for cognition: communication through coherence. Neuron 88, 220–235. doi: 10.1016/j.neuron.2015.09.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris, A. Z., and Gordon, J. A. (2015). Long-range neural synchrony in behavior. Annu. Rev. Neurosci. 38, 171–194. doi: 10.1146/annurev-neuro-071714-034111

PubMed Abstract | CrossRef Full Text | Google Scholar

Hurtado, J. M., Rubchinsky, L. L., and Sigvardt, K. A. (2004). Statistical method for detection of phase-locking episodes in neural oscillations. J. Neurophysiol. 91, 1883–1898. doi: 10.1152/jn.00853.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Izhikevich, E. M. (2007). Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. Cambridge, MA: MIT Press.

Lachaux, J. P., Rodriguez, E., Martinerie, J., and Varela, F. J. (1999). Measuring phase synchrony in brain signals. Hum. Brain Mapp. 8, 194–208. doi: 10.1002/(SICI)1097-0193(1999)8:4<194::AID-HBM4>3.0.CO;2-C

PubMed Abstract | CrossRef Full Text | Google Scholar

Le Van Quyen, M., and Bragin, A. (2007). Analysis of dynamic brain oscillations: methodological advances. Trends Neurosci. 30, 365–373. doi: 10.1016/j.tins.2007.05.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Oswal, A., Brown, P., and Litvak, V. (2013). Synchronized neural oscillations and the pathophysiology of Parkinson's disease. Curr. Opin. Neurol. 26, 662–670. doi: 10.1097/WCO.0000000000000034

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, C., Worth, R. M., and Rubchinsky, L. L. (2010). Fine temporal structure of beta oscillations synchronization in subthalamic nucleus in Parkinson's disease. J. Neurophysiol. 103, 2707–2716. doi: 10.1152/jn.00724.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, C., Worth, R. M., and Rubchinsky, L. L. (2011). Neural dynamics in parkinsonian brain: the boundary between synchronized and nonsynchronized dynamics. Phys. Rev. E 83:042091. doi: 10.1103/physreve.83.042901

PubMed Abstract | CrossRef Full Text | Google Scholar

Pikovsky, A., Rosenblum, M., and Kurths, J. (2001). Synchronization: A Universal Concept in Nonlinear Sciences. Cambridge: Cambridge University Press.

Google Scholar

Pittman-Polletta, B. R., Kocsis, B., Vijayan, S., Whittington, M. A., and Kopell, N. J. (2015). Brain rhythms connect impaired inhibition to altered cognition in schizophrenia. Biol. Psychiatry 77, 1020–1030. doi: 10.1016/j.biopsych.2015.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Ratnadurai-Giridharan, S., Zauber, S. E., Worth, R. M., Witt, T., Ahn, S., and Rubchinsky, L. L. (2016). Temporal patterning of neural synchrony in the basal ganglia in Parkinson's disease. Clin. Neurophysiol. 127, 1743–1745. doi: 10.1016/j.clinph.2015.09.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Rubchinsky, L. L., Ahn, S., and Park, C. (2014). Dynamics of desynchronized episodes in intermittent synchronization. Front. Comput. Phys. 2:38. doi: 10.3389/fphy.2014.00038

CrossRef Full Text | Google Scholar

Schnitzler, A., and Gross, J. (2005). Normal and pathological oscillatory communication in the brain. Nat. Rev. Neurosci. 6, 285–296. doi: 10.1038/nrn1650

PubMed Abstract | CrossRef Full Text | Google Scholar

Somers, D., and Kopell, N. (1993). Rapid synchronization through fast threshold modulation. Biol. Cybern. 68, 393–407. doi: 10.1007/BF00198772

PubMed Abstract | CrossRef Full Text | Google Scholar

Spellman, T. J., and Gordon, J. A. (2015). Synchrony in schizophrenia: a window into circuit-level pathophysiology. Curr. Opin. Neurobiol. 30, 17–23. doi: 10.1016/j.conb.2014.08.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Uhlhaas, P. J., and Singer, W. (2006). Neural synchrony in brain disorders: relevance for cognitive dysfunctions and pathophysiology. Neuron 52, 155–168. doi: 10.1016/j.neuron.2006.09.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Uhlhaas, P. J., and Singer, W. (2010). Abnormal neural oscillations and synchrony in schizophrenia. Nat. Rev. Neurosci. 11, 100–113. doi: 10.1038/nrn2774

PubMed Abstract | CrossRef Full Text | Google Scholar