Emergent gamma synchrony in all-to-all interneuronal networks

1. Introduction

Cortical gamma rhythms (30–80 Hz) are correlated with diverse brain functions such as memory formation (Singer and Gray, 1995), linguistic processing (Pulvermller et al., 1995), and associative learning (Ritz and Sejnowski, 1997; Miltner et al., 1999). Two commonly known mechanisms for gamma oscillations are Pyramidal-Interneuronal Network Gamma (PING) and Interneuron Network Gamma (ING).

Interneurons play critical roles for gamma rhythm generation in both mechanisms. However, in ING, interneurons are solely responsible for gamma frequency activity (Bartos et al., 2002). Additionally, there is evidence that interneurons show phase-locked activity with gamma oscillations (Bartos et al., 2007). In particular, perisomatic inhibition caused by fast-spiking parvalbumin (PV) immunoreactive basket cell interneurons (Lytton and Sejnowski, 1991) have been implicated in the generation of cortical and subcortical gamma oscillations (Bartos et al., 2007).

Previous theoretical works have shown that synchronization is possible in interneuronal networks (Wang and Rinzel, 1992; van Vreeswijk et al., 1994), however, this synchronization is not resilient to even slight heterogeneity or noise in the network (Wang and Buzsáki, 1996). In biological systems, inherent noise and heterogeneities in neuronal firing rates are likely to be present. Investigations by Kopell and Ermentrout (2004) have shown that gap junctions are one way by which interneuronal network synchronization can be robust against heterogeneities and noise in the network.

Evidence from literature suggests a strong correlation between Hebbian learning and Gamma rhythms (Miltner et al., 1999). Indeed, studies have established a link between glutamatergic synaptic plasticity and gamma oscillations (Whittington et al., 1997; Traub et al., 1998). More recently, synaptic plasticity of inhibitory Gamma-Aminobutyric acid (GABA) synapses was observed and characterized (Haas et al., 2006). This led us to question whether inhibitory synaptic plasticity could also play a role in generation of gamma oscillations.

In our earlier works, we demonstrated that inhibitory plastic synapses could improve in-phase synchronization in the presence of firing rate heterogeneity in a simple network of uni-directionally coupled interneurons (UCI) (Talathi et al., 2008). In the current paper, we investigate how the presence of synaptic plasticity can increase the range of firing rate heterogeneity for which mutually-coupled parvalbumin positive interneuronal (MCI) networks can synchronize at gamma frequencies. We begin with analyzing in-phase synchrony of a two interneuron MCI network (2-MCI) and prove that the analysis holds even for a much larger network (100-MCI). We demonstrate that the spike time response curve (STRC) of parvalbumin positive interneurons exhibit a strong second order component. We present a novel approach to account for the second order STRC component in the analysis of 1:1 in-phase synchronization of the 2-MCI parvalbumin positive interneuronal network.

The organization of the paper is as follows: the Methods Section presents the mathematical model for the network, the interneuron, and the synapse used in our studies. We analyze two types of heterogeneity in the interneuronal networks: (1) temporal heterogeneity; which corresponds to the heterogeneity in the intrinsic firing rates of the coupled interneurons, and (2) structural heterogeneity; which corresponds to the heterogeneity in the synaptic coupling strength. We introduce a few network measures in order to quantify the network state in terms of its connectivity and/or synaptic strengths. We then introduce a novel method for using STRC's with a strong second order component to estimate the spike times for a periodically firing neuron that receives periodic synaptic inputs. Finally, we introduce the empirical STDP rule (Haas et al., 2006; Talathi et al., 2008) used in our investigations of the emergence of synchronization in MCI networks.

In the Results Section, we systematically analyze the synchronization domain of a 2-MCI system in terms of temporal and structural heterogeneity. We show how heterogeneity and initial values of synaptic strengths determine the final synchronization state of the 2-MCI system. The 1:1 domain of synchronization identified by this analysis allows us to determine how STDP must evolve synaptic strengths in order to achieve 1:1 synchronization. We then use the STDP rule and STRC of the interneuron model to derive a nonlinear map for the evolution of the 2-MCI system to a state of in-phase synchronization. Finally, we demonstrate that the principal findings from our analysis of the 2-MCI network hold for larger MCI networks.

2. Materials and Methods

2.1. The Network Model

Each neuron in the network of all-to-all coupled inhibitory interneurons is modeled based on a single compartment model for parvalbumin positive inhibitory neuron developed by Wang and Buzsáki (1996) with a fast sodium channel, a delayed rectifier potassium channel and a leak channel. The dynamical equation for the model neuron is given by,

$\begin{array}{cc}\begin{array}{l}C\frac{d{V}_j(t)}{dt}=I_j^{DC}+g_{Na}{m}_{\infty }^{3}h(t)(E_{Na}-V_j(t))\hfill \\ +g_K{n}^4(t)(E_K-V_j(t))+g_L(E_L-V_j(t))\hfill \\ +{\displaystyle \sum _i{a}_{ij}}{g}_{ij}(t)s_i(t)(E_I-V_j(t))\hfill \end{array}\hfill & (1)\hfill \end{array}$

where C = 1 μ F/cm². V_j(t) is the membrane potential of the jth neuron (j = {1, 2, ⋯N}). The external current, $I_j^{DC}$ is given as follows:

$\begin{array}{cc}{I}_j^{DC}=I_{ref}+\left(j-\frac{N+1}{2}\right)\left(\frac{H\times I_{ref}/100}{N-1}\right)& (2)\end{array}$

where, the parameter H controls the degree of heterogeneity in the intrinsic firing activity of the interneurons in the network. When H > 0, $I_j^{DC}<I_{j+1}^{DC}$. Since neuron j's firing period is determined by $I_j^{DC}$, this means that neuron j fires slower than neuron j + 1. The actual value of H determines the amount by which neighboring neurons differ in their I^DC-values (and correspondingly, firing periods), as given in Equation (2). Hence for a larger H, there is a larger range of firing periods across the neurons in the network, corresponding to a larger temporal heterogeneity. The parameter I_ref = 1 μA/cm² in Equation (2) is set such that the mean intrinsic frequency of firing for the neurons in the network is ≈ 60 Hz (Wang and Buzsáki, 1996). E_r (r = Na, K, L) are reversal potentials of the sodium and potassium ion channels and the leak channel, respectively. E_I, is the reversal potential of the fast GABAergic inhibitory synapse. g_r (r = Na, K, L) represent the conductance of sodium, potassium, and the leak channel, respectively. The steady state activation for sodium current, m_∞ = α_m∕(α_m + β_m). The inactivation variable for sodium channel h(t) and the activation variable for potassium current n(t) satisfy the following first order kinetic equation: $\frac{dX(t)}{dt}=\varphi ({\alpha }_X(V_j(t))(1-X(t))-{\beta }_X(V_j(t))X(t))$, where X(t) = h(t), n(t) with ϕ = 5. The functions α_X and β_X are given by:

$\begin{array}{cc}\begin{array}{l}\begin{array}{cc}{\alpha }_m(x)=\frac{-0.1(x+35)}{e^{-0.1(x+35)}-1}\hfill & {\beta }_m(x)=4e^{-(x+60)/18}\hfill \end{array}\hfill \\ \begin{array}{cc}{\alpha }_h(x)=0.07e^{-(x+58)/20}\hfill & {\beta }_h(x)=\frac{1}{e^{-0.1(x+28)}+1}\hfill \end{array}\hfill \\ \begin{array}{cc}{\alpha }_n(x)=\frac{-0.01(x+34)}{e^{-0.1(x+34)}-1}\hfill & {\beta }_n(x)=\frac{0.125}{e^{-(x+44)/88}}\hfill \end{array}\hfill \end{array}\hfill & (3)\hfill \end{array}$

The strength of inhibitory synaptic coupling from the pre-synaptic neuron i onto the post-synaptic target neuron is modeled via the variable g_ij(t). In the presence of spike timing dependent plasticity; the coupling strength evolves depending on the plasticity rule as a function of the relative timing of spiking activity between the neurons i and j. In absence of synaptic plasticity the synaptic coupling strength is set as:

$\begin{array}{cc}{g}_{ij}=\frac{g_0}{N}\left(1+\frac{\eta }{100}\text{sgn}(i-j)\right)& (4)\end{array}$

where |η| ≤ 100. For all the results presented unless otherwise stated, g₀ = 0.1 mS/cm². The parameter η regulates the structural heterogeneity in the network such that for η ≠ 0, there is asymmetry in the synaptic coupling strengths of the two mutually coupled neurons. By design, the asymmetry is such that for η > 0, the strength of synapse from neuron with lower index is smaller than that of the synapse from neuron with higher index label. The variable a_ij defines the topology of the inhibitory neuronal network such that, for an all-to-all coupled network, we have a_ij = 1 ∀i ≠ j and a_ij = 0 otherwise for all i, j = {1, 2, ⋯N}. The variable s_i(t) models the fraction of bound neurotransmitters onto the post-synaptic neuron GABAergic receptors resulting from a spike elicited by the pre-synaptic neuron i, and satisfies the following first order kinetic equation, Abarbanel et al. (2003) and Talathi et al. (2008),

$\begin{array}{cc}{\dot{s}}_i(t)=\frac{S_0(V_i(t))-s_i(t)}{\widehat{\tau }(S_I-S_0(V_i(t))}& (5)\end{array}$

The kinetic equation for s_i(t) involves two time constants, ${\tau }_R=\widehat{\tau }(S_I-1)=0.1$ ms, the docking time for the neurotransmitter binding and ${\tau }_D=\widehat{\tau }{S}_I=10$ ms, the undocking time constant for the neurotransmitter binding. Finally, S₀(x) is the sigmoidal function given by, S₀(x) = 0.5(1 + tanh(120(x − 0.1))).

In order to quantify the influence of synaptic plasticity on the structure of synaptic connectivity in the network, we define the following network measures:

(a) Link imbalance: The difference in synaptic strengths of the two coupled neurons,

$\begin{array}{cc}{L}_{ij}=-L_{ji}=a_{ij}{g}_{ij}-a_{ji}{g}_{ji}& (6)\end{array}$

(b) Neuronal strength: The sum of synaptic strength of all outgoing synapses from the neuron i,

$\begin{array}{cc}{G}_i={\displaystyle \sum _j{a}_{ij}}{g}_{ij}& (7)\end{array}$

For all the networks considered here, the neurons in the network are labeled in ascending order of input DC current such that label 0 is assigned to the neuron with lowest intrinsic firing rate while label N − 1 is assigned to neuron with the highest intrinsic firing rate.

2.2. Spike Time Response Curve

As a measure of the influence of synaptic input on the firing times of a given neuron i, we define the spike time response curves (STRC's) ${\Phi }_{ij}({\tau }_R,{\tau }_D,E_R,T_{i0},g,\delta t)=\frac{T_{ij}-T_{i0}}{T_{i0}}$ (Acker et al., 2004; Oprisan et al., 2004; Talathi et al., 2010), where T_i0 is the intrinsic period of spiking, and T_ij represents the length of the jth spiking cycle from the cycle j = 1 in which the neuron receives synaptic input at time 0 < δt <T_i0. A synaptic perturbation can either advance or delay the occurrence of an impending spike depending on the bifurcation character of the neuron (Ermentrout, 1996) and the type of synaptic perturbation (excitatory vs. inhibitory). The number of firing cycles j, for which the effect of synaptic perturbation lasts, depends on the synaptic parameters: the synaptic rise time τ_R, the synaptic decay time τ_D, the reversal potential of the synapse E_R, and the synaptic conductance g. For fast hyperpolarizing inhibitory synapse (E_R = −75 mV), as considered in this study, a synaptic input delays the time of occurrence of subsequent spike such that Φ_i1 ≥ 0 ∀δt > 0, Φ_i2 > 0 for δt → T_i0, and Φ_ij = 0 ∀j > 2 (Talathi et al., 2010). We note that in all further calculations, unless otherwise mentioned, we suppress the dependence of Φ on τ_R, τ_D, E_R, and the intrinsic firing period of neuron T₀ and define Φ(τ_R, τ_D, E_R, T₀, g, δt) ≡ Φ(g, δt).

We begin by estimating STRC's using the direct method as follows: consider neuron i, firing regularly with period T_i0. The neuron is perturbed through an inhibitory synapse at time δt after the neuron has fired a spike at reference time zero (see Figure 1A). The spiking time for a neuron is considered to be the time when its membrane voltage V, crosses a threshold (set to 0 mV in all the calculations presented here). As a result of this perturbation, the neuron fires the next spike at time t₁, representing the first cycle after perturbation of length T_i1 ≠ T_i0. Depending on the properties of the synapse, i.e., g_s, τ_R, τ_D, and E_R; the length of subsequent cycles might change. In Figure 1B, we plot the first two STRC's Φ_i1 (solid black line) and Φ_i2 (solid magenta line), estimated by using the procedure described above for δt ∈ [0, T_i0). We note the presence of a non-zero second order STRC component for the choice of the parameters of the neuron.

FIGURE 1

Figure 1. STRC Calculations. (A) Schematic for determining the STRC using the standard method. (B) The neuron firing intrinsically at 60 Hz, receives inhibitory synaptic perturbation with synaptic parameters: E_R = −75 mV, τ_R = 0.1 ms, τ_D = 5 ms, and g = 0.1 mS/cm². The original STRC is shown in solid colors while the STRC with the included phase correction term (δt*) is represented by the dotted lines. (C) Schematic of an uni-directionally coupled two-interneuron network. (D) For H = 30 the neurons exhibit phase-locked 1:1 synchronization. The spike from Neuron 0 (black voltage trace) synaptically perturbs neuron 1 (magenta voltage trace) at a fixed δ_t during every firing cycle of neuron 1. The solid blue lines indicate the spike-time shift calculated using the STRC that incorporates the correction factor. The dotted blue lines indicate the spike-time shift calculated using the original STRC's without any correction term.

In the direct method for computing STRC's, we use a single perturbation to observe the changes in spike times of a periodically firing neuron. These changes in spike times can be used to predict future spike times. If Φ_i2 ≈ 0, then a single perturbation arriving at δt, following the last spike at reference time t₀, can be assumed to produce a spike time shift such that the next spike occurs at time t₁ = t₀ + T_i0(1 + Φ_i∞(g, δt)), where Φ_i∞(g, δt) = Φ_i1(g, δt) + Φ_i2(g, δt) (Oprisan and Canavier, 2001). When a significant higher order STRC component (Φ_i2) is present, Φ_i∞(g, δt) can no longer be expressed as a linear sum of the first and second order STRC components (Talathi et al., 2010). With a significant second order STRC component, the neuron requires at least two firing cycles before it resets to its original firing period. In such conditions, the assumption of linear second order STRC correction to the estimated future spike time is no longer valid (Talathi et al., 2010).

We demonstrate this issue using a 2-UCI network as shown in Figure 1C. Here, neuron 1, which fires periodically, is also periodically perturbed by neuron 0 via an inhibitory synapse. We consider a specific case of the 2-UCI network where firing rate heterogeneity H = 30 (Equation 2) and phase locked 1:1 synchrony exists. The voltage traces of neuron 0 and 1 are shown as black and magenta solid curves, respectively. We see in Figure 1D that when the shift in spike times are predicted using T_i1 = T_i0(1 + Φ_i∞(g, δt)), there is a noticeable error between predicted spike times (blue dotted stems) and the actual spike times (black solid curve). We next sought to empirically quantify the error between predicted and actual spike times as a function of δt.

Previously, Talathi et al. (2008) derived a STRC based theoretical map for when a 2-UCI network (see Figure 1C) exhibits steady state 1:1 phase locked synchrony. The map is expressed as:

$\begin{array}{cc}{T}_{00}-T_{10}=T_{10}{\Phi }_{1\infty }(g_{01},\delta t_s)& (8)\end{array}$

where, T₀₀ and T₁₀ represent the firing periods of neurons 0 and 1, respectively, Φ_1∞(g₀₁, δt_s) is the STRC calculated for the steady-state perturbation time δt_s, which is effectively the spike time difference between the two neurons. The spike time difference could also be measured directly by observing simulation results of the 2-UCI network in a state of 1:1 synchrony, which we denote as δt*. When the STRC's second order component is weak, the theoretical predictions match closely with the simulations and we get $\delta t^{*}\approx \delta t_s$. However, when the STRC's second order component is strong, there is a noticeable deviation between the map's δt_s and simulation's δt*. We term this error as $\overline{\delta t}=\delta t^{*}-\delta t_s$.

In order to characterize $\overline{\delta t}$ as a function of δt_s, we need to observe the 2-UCI network in 1:1 synchrony across a range for δt* and δt_s. Since, δt_s and δt* are spike time differences of the theoretical map solution and simulations, respectively for a particular configuration of the 2-UCI network in 1:1 synchrony, they cannot be arbitrarily predetermined. Instead, by varying the temporal heterogeneity parameter H (25 ≤ H ≤ 45), we can cause the 2-UCI network to exhibit a range of 1:1 phase-locked synchrony states. Each of these 1:1 synchrony states allows us to measure the true spike time differences, δt*, from simulations and the theoretically estimated spike time difference, δt_s, between the two neurons. By noting δt_s, δt*, and the corresponding error $\overline{\delta t}$ for every value of H that produces 1:1 in-phase synchronization in the 2-UCI network, we were able to plot the error ($\overline{\delta t}$) as a function of δt_s (see Figure 2A). We noticed that the relationship between $\overline{\delta t}$ vs. δt_s is approximately linear and can be approximated via a linear fit as $\overline{\delta t}=-0.035064\times \delta t_s+2.2955$. This fitted line allowed us to calculate the approximate correction factor that has to be added to a given δt_s in order for the STRC based solution to match simulation results (see Figure 2B). By generalizing δt_s to any perturbation δt, the correction factor can be expressed as $\overline{\delta t}=-0.035064\times \delta t+2.2955$. The correction factor ($\overline{\delta t}$) then lets us continue expressing the STRC as the linear sum of its first and second order components:

$\begin{array}{cc}{\Phi }_{i\infty }(g,\delta t^{\ast })={\Phi }_{i1}(g,\delta t^{\ast })+{\Phi }_{i2}(g,\delta t^{\ast })& (9)\end{array}$

$\begin{array}{cc}\delta t^{\ast }=\delta t+\overline{\delta }t=0.964936\ast \delta t+2.2955& (10)\end{array}$

FIGURE 2

Figure 2. 2-UCI Simulation and STRC Based Spike Time Calculations. (A) Shows the variation of simulation spike time difference δt* and STRC based map predicted δt_s as H is varied. The negative values for δt_s is due to the inaccuracy of the standard STRC method arising from second order STRC components. (B) Shows the difference between actual and estimated spike time differences in the 2-UCI network. A linear fit approximately describes the relationship between $\overline{\delta t}$ and δt_s as: $\overline{\delta t}=-0.{035064}^{*}\delta t_s+2.2955$.

We will assume that all further references to STRC's in this paper incorporate the correction term δt* that is calculated from the given δt.

2.3. Spike Timing Dependent Plasticity Rule

We will consider the spike timing dependent plasticity (STDP) rule for GABAergic synapses, identified at GABAergic synapses in the entorhinal cortex (Haas et al., 2006), which is similar to the well-known STDP rule for excitatory synapses (Bi and Poo, 1998). The STDP rule is asymmetric with long-term potentiation for positive time intervals, Δt, between the timing of post-synaptic neuron firing and the pre-synaptic neuron firing and long-term depression for Δt < 0. We assume Δt to be positive if the post-synaptic spike occurs after a pre-synaptic spike. The empirical function defining the STDP rule is given as:

$\begin{array}{cc}\begin{array}{l}\Delta g(\Delta t)=g+\tilde{g}(\Delta t)\text{if }\Delta t>0\hfill \\ =g-\tilde{g}(\Delta t)\text{if }\Delta t<0\hfill \end{array}\hfill & (11)\hfill \end{array}$

where $\tilde{g}(\Delta t)=\frac{1}{{\beta }^{\beta }{e}^{-\beta }}{\alpha }^{\beta }\left|\Delta t\right|\Delta t^{\beta -1}{e}^{-\alpha \left|\Delta t\right|}$ and g₊ ≥ g₋. The parameters are β = 10 and α = 0.94, which provide a temporal window of ±20 ms over which the efficacy of the synaptic plasticity is non-zero (Haas et al., 2006). Figure 3 presents the learning rule we use for STDP. In all the results presented here, unless otherwise mentioned, we set g₊ = g₋ = 0.01 and we assume an additive (linear) update rule for the modification of the synaptic strengths in the network (Caporale and Dan, 2008). Furthermore, we implement the STDP rule through all-to-all interactions, whereby after a spike occurs in a given post-synaptic neuron all the synapses from pre-synaptic neurons with immediately preceding spikes are potentiated whereas all the synapses from pre-synaptic neurons with spikes occurring after the post-synaptic spike are depressed (Bi and Wang, 2002; Froemke and Dan, 2002; Froemke et al., 2006).

FIGURE 3

Figure 3. Schematic of the spike timing dependent plasticity rule. The parameters are: α = 0.94, β = 10, g₊ = g₋ = 0.01 mS/cm².

2.4. Synchronization Metric

We will measure the degree of synchronization S in the N-MCI networks by quantifying the amount of fluctuations in the network activity (Hansel and Sompolinsky, 1992; Kudela et al., 2003). Given the values for the membrane potentials V_i(t − kt_s) where i = (0, 1, ⋯N − 1), k = (0, 1, ⋯M − 1), t_s is the sampling time interval and M quantifies the time horizon T = Mt_s, we define:

$\begin{array}{cc}S(N)=\frac{N{\sigma }_V}{{\displaystyle {\sum }_{i=0}^{N-1}{\sigma }_{V_i}}}& (12)\end{array}$

where ${\sigma }_V^2$ represents the variance in the mean field activity of the network and ${\sigma }_{V_i}^2$ represents the variance in the membrane potential of ith neuron in the network measured over a large time interval T. If we define $\langle X{\rangle }_t=\frac{1}{T}{\int }_{t-T}^{t}X(t)dt$, then σ_V and σ_Vi are given as:

$\begin{array}{l}{\sigma }_{V_i}^2={\langle V_i^2\rangle }_t-{\langle V_i\rangle }_t^2\hfill \\ {\sigma }_V^2={\langle {\left(\frac{1}{N}{\displaystyle \sum _i{V}_i}\right)}^2\rangle }_t-{\langle \frac{1}{N}{\displaystyle \sum _i{V}_i}\rangle }_t^2\hfill \end{array}$

The properties of S are well-characterized in literature (Golomb and Rinzel, 1993). S essentially measures the average fluctuations of the mean membrane voltage of the network over a long time duration. It is normalized by the sum of individual neuron membrane fluctuations such that S scales between 0 and 1. In particular, for a fully synchronized network such that V_i(t) = V(t) ∀i ∈ [0, N − 1], S = 1, whereas 0 ≤ S < 1 in all other situations.

3. Results

3.1. Synchronization Manifold of the 2-MCI Network

In this section, we investigate the structure of the synchronization manifold of the 2-MCI network as a function of the network heterogeneity parameters: {H, η}. We emphasize that all of our results critically depend on the structural heterogeneity parameter, η (Equation 4) and on the temporal heterogeneity parameter H (Equation 2). Now, let t_{0, i} (i = 1, 2, ⋯ ) and t_{1, j} (j = 1, 2, ⋯ ) represent the spike times for the two coupled neurons 0 and 1 in the 2-MCI network (Figure 4A). We define the mean period $\langle T_X\rangle =\underset{n\to \infty }{lim}\frac{t_{X,n}}{n}$ (X = 0, 1) and the time lag δ_n = t_1,n − t_0,m between the nearest spike times of the two neurons. The index m corresponds to the mth spike from neuron 0 such that t_1,n−1 ≤ t_{0, m} ≤ t_1,n. For a given value of H and η, we numerically solve Equation (1) in order to estimate 〈T₀〉 and 〈T₁〉. If there exists m, n ∈ ℕ and m ≠ 0 and n ≠ 0 such that the ratio:

$\begin{array}{cc}R=\frac{\langle T_0\rangle }{\langle T_1\rangle }=\frac{m}{n}& (13)\end{array}$

we label the point {H, η} on the two dimensional hyper-plane spanned by H and η as the m:n synchronization point. Hence R describes the order of a synchronized pair of neurons as the ratio of their firing periods. The results of these calculations for 0 ≤ H < 50 and −50 ≤ η < 50 are summarized in Figure 4A. For example, the black dots in Figure 4A span the region of 1:1 synchronization. In this situation δ_n approached a fixed value δ_s ≤ 〈T_X〉 representing a phase-locked state (see Figures 4B–D). The regions of higher order synchronization are depicted in various colors as follows: R = 6:5 (in orange); R = 4: 3 (in green); R = 3: 2 (in blue); R = 5: 3 (in magenta); and R = 2: 1 (in cyan). An important feature of the structure of the synchronization manifold is that the range of values of dynamical heterogeneity H over which the 2-MCI network can sustain 1:1 synchronization is significantly enhanced as the value of structural heterogeneity η decreases below zero. In other words, the effect of increasing H, resulting in an increase in the mismatch between the intrinsic firing rates of the two coupled neurons, can be compensated by increasing the strength (S) and simultaneously decreasing the sensitivity (K) of the neuron with low intrinsic firing rate. Yet another feature of the synchronization manifold is that the minimum absolute value of structural heterogeneity |η| required to sustain 1:1 synchronization in the 2-MCI network increases monotonically with increasing heterogeneity H. Furthermore, islands of higher order synchronization states emerge as H increases for a given fixed value of η. We will see in Section 3.2, this property of the synchronization manifold will significantly influence the evolution of the 2-MCI network under the influence of STDP.

FIGURE 4

Figure 4. Synchronization manifold of the 2-MCI network. (A) Synchronization manifold of the 2-MCI network (shown in inset). The Figure shows the domain of m:n synchronization for the 2-MCI network in the two dimensional plane spanned by the network heterogeneity parameter H, and the network imbalance parameter η. In (B–D), we plot the time lag δ between the firing times of the two neurons in the 2-MCI network as function of the network temporal heterogeneity parameter H for three different cases of the network temporal heterogeneity parameter η = {−20, 0, 20}, respectively. In each Figure, we also plot the ratio of mean firing periods of the two neurons R (magenta dots).

In Figures 4B–D, we plot δ_n as function of H for the cases η = {−20, 0, 20}, respectively. We also plot R (magenta dots in Figures 4B–D, y-axis on the right) as a function of H. We see from Figure 4C that for a perfectly homogeneous 2-MCI network, i.e., H = η = 0, the two coupled neurons exhibit 1:1 in-phase synchronization. In the absence of structural heterogeneity (η = 0), the network can sustain 1:1 synchronization for 0 ≤ H < 8. However, as can be seen from Figure 4A, for η = −20, the range of values for H over which the network can sustain 1:1 synchronization is enhanced to 1 < H < 15. When η < −30, we note that the range of H-values for 1:1 synchronization does not include H = 0, which suggests that the structural heterogeneity is such that the network tends to exhibit master-slave dynamical structure wherein the influence of one neuron (neuron 0 in this case) on to the target neuron (neuron 1 in this case) is much stronger than in the vice versa case. In the case when η = 20, the bias toward master-slave dynamics is significantly stronger with neuron 1 being the driver neuron and neuron 0 being the slave neuron. This is because, not only is neuron 1 firing at an intrinsically greater rate than neuron 0, but also the synaptic drive onto neuron 0 is stronger, which further slows the firing rate of neuron 0. The network is therefore unable to exhibit 1:1 synchronization. At the same time, the domain of higher order synchronization such as 2:1 synchronization is significantly enhanced.

In summary, the key conclusion from this analysis is that the 2-MCI network can sustain 1:1 synchronization for a wide range of heterogeneity in the intrinsic firing rates of the coupled neurons for an appropriate choice of structural heterogeneity in the network, i.e., by increasing strength of neuron 0 relative to the strength of neuron 1 and decreasing the sensitivity of neuron 0 relative to that of neuron 1, to trigger an effective structural imbalance η < 0.

3.2. STDP Induced Emergent Synchronization in 2-MCI Network

We will now investigate the effect of STDP on the dynamics of a heterogeneous 2-MCI network. For all the results presented in this section, we will begin with a structurally homogeneous network i.e., initial η = η₀ = 0 and investigate the influence of STDP on the network dynamics for various values for H. For a given value of 0 ≤ H < 50, we numerically solve Equation (1) beginning with a random set of initial conditions V_i (i= 0, 1) in the range [−70, −50] mV. We perform simulations (each simulation run is for 5 s) for each of the two cases: (1) no STDP learning, with static synaptic conductance values g₀₁ = g₁₀ = g₀∕2 and (2) with STDP learning, following the additive rule given in Equation (11). The initial synaptic conductance values are set to be the same as for the static case. For each simulation run, after the network has evolved into an asymptotic state (asynchronous firing or the state of m:n synchronization) we estimate the ratio R (see Equation 13). The results of these calculations are summarized in Figure 5, where we plot R as function of H and the color code represents the fraction P of the simulations generating a given value of R. For example, from Figure 5A, we see that the 2-MCI network exhibits a stable state of 1:1 synchronization for 0 ≤ H < 9. The network also exhibits narrow bands of stable higher order synchronization for larger heterogeneity values (as seen in Figure 4A). It is clear from Figure 5A that many of these states of stable synchrony are global attractors for the network, i.e., the state of stable synchronization is achieved independent of the initial conditions (P = 1) on the membrane potentials of the neurons in the network.

FIGURE 5

Figure 5. 2-MCI Synchronization Ratio probability for η₀ = 0. (A) Shows the probability of the 2-MCI network to reach different orders of synchronization without STDP. (B) Shows the probability for different orders of synchronization with STDP. The range of H for which the network can synchronize is significantly enhanced for the 1:1 domain and higher orders.

3.2.1. STDP Induced 1:1 Synchronization

In the presence of STDP, the synaptic coupling strengths evolve such that, the range of heterogeneity values over which the network can sustain stable states of 1:1 synchronization (P = 1) is enhanced to include the range 0 ≤ H < 24 (see Figure 5B). In Figure 6A, we plot the mean value (and the standard error) of synchronization metric S as function of H. We note that S = 1 for the range of H-values over which the network can sustain a globally stable attractor state of 1:1 synchronization, representative of the case when the network is fully synchronized (i.e., 1:1 in-phase synchronization). As noted in Section 3.1, in absence of STDP, only the homogeneous 2-MCI network can exhibit fully synchronized state (also see the black curve in Figure 6A). Finally, as can be seen from Figure 5A, while the likelihood for the 2-MCI network to exhibit 1:1 synchronization is zero (P = 0) for H > 8, with STDP learning (Figure 5B) the likelihood for the 2-MCI network to sustain 1:1 synchronization, P ≠ 0, for the range of 0 ≤ H < 27, albeit the likelihood for the network to sustain 1:1 synchronization decreases with increasing H.

FIGURE 6

Figure 6. 2-MCI Synchrony metric. (A) A measure of the synchrony metric when STDP is enabled (black) and without STDP (magenta) for the range 0 ≤ H ≤ 50. (B) Shows the evolution of η for an example initial case of {η₀ = 0, H = 10}. In the inset we show the corresponding voltage traces of the neurons after reaching in-phase synchronization.

In order to understand how the structural properties of the 2-MCI network change as the network evolves to a state of 1:1 in-phase synchronization under the influence of STDP, in Figure 6B, we plot the time evolution of η for the case H = 10. The inset of the Figure shows a sample time-trace of the membrane potential of the two neurons after the network has evolved into the asymptotically stable state of 1:1 in-phase synchronization. The network is initially structurally symmetric, i.e., η₀ = 0. The STDP learning begins at time t = 200 ms, immediately after which, the network imbalance begins to evolve ultimately reaching a steady state value η ≈ −40.

Starting from η₀ = 0, the network evolves vertically downwards in terms of the two-dimensional plot of the synchronization manifold in Figure 4A, such that the synaptic strength of the slower firing neuron is enhanced while at the same time its sensitivity is decreased. As noted earlier (Section 3.1), this is the type of structural heterogeneity required for the 2-MCI network to sustain 1:1 synchronization for H > 8. This evolution of the 2-MCI network to a preferred state of 1:1 in-phase synchronization is critically dependent on the form of the learning rule that corresponds to Δg > 0(< 0) for Δt > 0(< 0), respectively. This form of the learning rule implies that for a given synapse, every time the post-synaptic neuron fires, the strength of synapse will increase where as every time the pre-synaptic neuron fires, the synaptic strength will decrease. For the configuration of the 2-MCI network with H = 20, neuron 1 has a greater intrinsic firing rate than neuron 0 such that in absence of STDP, on average 〈T₁〉 < 〈T₀〉. Thus, after STDP is active, the synapse g₀₁ (with neuron 1 being the post-synaptic neuron for the synapse) increases more often than it decreases and consequently η fluctuates with an increased bias toward η < 0 until the network evolves into the domain of 1:1 synchronization (see Figure 4A) at which point 〈T₁〉 = 〈T₀〉. Since Δg ≠ 0 for Δt ≠ 0, the network continues to evolve toward the fixed point δ = 0 corresponding to the stage where η monotonically decreases (Figure 6B) until the network reaches a stable state of 1:1 in-phase synchronization. Hence in the case of 1:1 synchrony, STDP evolves the structural heterogeneity (η) such that it compensates for any temporal heterogeneity (H > 0).

3.2.2. STDP Induced m:n Synchronization

We can see from Figure 5B that for H > 24, under the influence of STDP, the state of 1:1 synchronization is not a global attractor of the 2-MCI network, i.e., P < 1. The network can evolve to states of higher order synchronization depending on the initial conditions on the membrane potentials of the neurons. It is also clear from Figure 5B that the likelihood (i.e., the fraction of simulation runs) for the network to evolve to a state of higher order synchronization increases with increasing H. For example, for H > 45, the probability for the network to evolve to a state of stable 2:1 is 1.

The synchronization manifold of the static 2-MCI network (Figure 4A) provides an insight into the reason for 1:1 synchronization not being the global attractor of the 2-MCI network under the influence of STDP learning. As discussed in Section 3.2.1, the STDP rule modulates the synaptic strengths such that the 2-MCI network evolves from a structurally balanced state of η₀ = 0 toward a structurally imbalanced state with η < 0 in order to compensate for H ≠ 0. However, for H ≥ 24, the synchronization manifold of the 2-MCI network exhibits multiple domains of higher order synchronization states for η < 0 and there is a likelihood for the network to enter these domains of higher order synchronization before the network can reach the state of 1:1 synchronization. As H increases beyond H = 18, the number of higher order synchronization states increase and at the same time the domain of 1:1 synchronization moves further away from the balanced state of η₀ = 0 and correspondingly the likelihood of the network to evolve to 1:1 synchronization decreases, as is evidenced from Figure 5B.

3.2.3. Analysis of Emergent 1:1 In-phase Synchrony in the 2-MCI Network

In this section, we analyze the stability of the emergent 1:1 in-phase synchronization in the 2-MCI network under the influence of STDP by deriving a nonlinear map for the evolution of the time lag δ between successive spike times of the two neurons in the 2-MCI network using the framework of STRC's. We consider the specific case of the network configuration with parameters: H = 10 and η₀ = 0. Following from the results presented in Section 3.1, for the choice of network heterogeneity and imbalance parameters, in the absence of STDP, the two neurons in the 2-MCI network will fire asynchronously, with mean firing rate of neuron 1 greater than that of neuron 0. Invoking STDP in this situation will, on average, cause the strength S₀ of neuron 0 to increase, while at the same time it will cause the sensitivity K₀ of neuron 0 to decrease. This is because, more often than not, the fast firing neuron 1 will emit more than one spike during each period of spiking from neuron 0. Each firing of neuron 1 (the post-synaptic neuron to synapse from neuron 0 to neuron 1), in turn, increases the strength of synapse g₀₁ through the STDP rule. By the same token, the strength of synapse g₁₀ will decrease creating an effective network imbalance η < 0. Thus, starting from a structurally balanced network I_net = 0, the network will evolve vertically downwards in the two dimensional η − H plane. As seen from Figure 6, for the specific case of H = 10 beginning with η₀ = 0, the network will eventually evolve into the domain of 1:1 synchronization.

In order to understand how the synaptic strengths evolve within the domain of 1:1 synchronization to the final stable state of 1:1 in-phase synchronization, we use the mathematical framework of STRC's to derive a nonlinear map for the evolution of the time-lag δ and the synaptic strengths g_ij (i, j = {0, 1}). We assume that once the network enters the phase locked state of 1:1 synchronization with 〈T₀〉 = 〈T₁〉, the phase locked state remains quasi-static as the synaptic strength continues to evolve toward the asymptotic state of 1:1 in-phase synchronization. In Figure 7, we present a snapshot of sequence of spikes from the two neurons after the neurons have entered the domain of 1:1 synchronization. Following from Figure 7, if t_0,n and t_1,n represents the timing of nth spike originating from neuron 0 and neuron 1 in the 2-MCI network, respectively, then the subsequent spike from the two neurons will emerge at times: t_0,n+1 and t_1,n+1 such that:

$\begin{array}{cc}\begin{array}{l}{t}_{0,n+1}=t_{0,n}+T_{00}\left[1+{\Phi }_{0\infty }\left(g_{10,n},(t_{1,n}-t_{0,n})\right)\right]\hfill \\ t_{1,n+1}=t_{1,n}+T_{10}\left[1+{\Phi }_{1\infty }\left(g_{01,n},(t_{0,n+1}-t_{1,n})\right)\right]\hfill \\ g_{01,n+1}={\tilde{g}}_{01,n}+\Delta g(t_{1,n+1}-t_{0,n+1})\hfill \\ g_{10,n+1}={\tilde{g}}_{10,n}+\Delta g(t_{0,n+1}-t_{1,n})\hfill \end{array}\hfill & (14)\hfill \end{array}$

where ${\tilde{g}}_{01,n}=g_{01,n}+\Delta g(t_{1,n}-t_{0,n+1})$ and ${\tilde{g}}_{10,n}=g_{10,n}+\Delta g(t_{0,n}-t_{1,n})$, respectively. In writing Equation (14), we assume that g_{ij, n} is the conductance of synapse from the pre-synaptic neuron i to the post-synaptic neuron j, immediately after nth spike is emitted from the post-synaptic neuron. We further assume that STDP induced change in synaptic conductance is instantaneous. The nonlinear map for the evolution of the time-lag δ_n = t_1,n − t_0,m (index m corresponds to the mth spike from neuron 0 such that t_{1, n−1} ≤ t_{0, m} ≤ t_{1, n}) and the synaptic conductances g_{ij, n} can be obtained from Equation (14) as:

$\begin{array}{cc}\begin{array}{l}{\alpha }_n=T_{00}\left[1+{\Phi }_{0\infty }(g_{10,n},{\tilde{\delta }}_n)\right]-{\tilde{\delta }}_n\hfill \\ {\delta }_{n+1}=T_{10}\left[1+{\Phi }_{1\infty }(g_{01,n},{\tilde{\alpha }}_n)\right]-{\tilde{\alpha }}_n\hfill \\ g_{01,n+1}=g_{01,n}+\Delta g(-{\alpha }_n)+\Delta g({\delta }_{n+1})\hfill \\ g_{10,n+1}=g_{10,n}+\Delta g(-{\delta }_n)+\Delta g({\alpha }_n)\hfill \end{array}& (15)\hfill \end{array}$

where $\tilde{\delta }=\delta mod{T}_{00}$ and $\tilde{\alpha }=\alpha mod{T}_{10}$.

FIGURE 7

Figure 7. 1:1 Phase locked synchrony in a 2-MCI network. This figure illustrates the evolution of consecutive spike times {δ, α} during phase locked 1:1 synchronization. The corresponding changes in {g₀₁, g₁₀} due to the STDP rule are also shown. We use this concept of 1:1 phase locked spike-time and STDP evolution to derive the map defined in Equation (15).

To determine whether Equation (15) can predict 1:1 in phase synchronization for the 2-MCI network under the influence of STDP learning, we apply the discrete map for the specific cases of the 2-MCI network with heterogeneity H = 10 and H = 20 with initial conditions g_{01, 0} = g_{10, 0} = g₀∕2 corresponding to the case η₀ = 0 and δ₀ = 0.

The results are presented in Figure 8. For the case H = 10, we see from Figure 8 that η evolves to a steady state value of η = −44, with δ_s = 18.8, α_s = 0.1 ms, g₀₁ = 0.073 mS/cm² and g₁₀ = 0.027 mS/cm². The mean period of the synchronized network (described in Section 3.1) is 〈T₀〉 = 〈T₁〉 = 18.9 ms. For the special case of 1:1 synchrony, we will refer to the network period as:

$\begin{array}{cc}\langle T_{0/1}\rangle =\langle T_0\rangle =\langle T_1\rangle & (16)\end{array}$

Upon simplifying (Equation 15) for the 1:1 steady-state, one can derive a straightforward relationship between δ_s and 〈T_0∕1〉 as: 〈T_0∕1〉 = δ_s + α_s. We see that when α_s ≈ 0 then 〈T_0∕1〉 = δ_s which is the in-phase 1:1 solution. Though the discrete map correctly predicts the existence of a steady state fixed point solution, a minor discrepancy in the discrete map solution for steady state value of η exists. There is a difference in η of ≈4 between the map's solution and to the solution obtained by numerically solving Equation (1) under the influence of STDP. This discrepancy stems from the linear approximation of the correction factor $\overline{\delta t}$ as a function of δt, as described in Section 2.2 regarding STRC's.

FIGURE 8

Figure 8. Map evolution of η. We validate that under STDP, the map evolves to stationary values of η that allow in-phase synchronization. For H = 10 the map predicts η = −44 (as compared to η = −40 for 2-MCI simulations). For H = 20 the map predicts η = −82 (as compared to η = −80 for 2-MCI simulations).

For the case H = 20, the steady state solution for the discrete map evolves to an in-phase solution at η = −82. As H approaches 27 and η = 0, the map fails to evolve to a steady state solution as expected. The map cannot accurately predict higher order synchronization, as Equation (15) is specifically derived for the 1:1 phase locked synchronization domain. From these findings we conclude that the discrete map for the temporal evolution of the spike time lag δ and the structural heterogeneity η is able to correctly predict the existence of a stable fixed point solution for 1:1 synchronization of the 2-MCI network under the influence of STDP learning.

We next determine if the 1:1 in phase is stable. First, we performed linear stability analysis of the discrete map in Equation (15) for the fixed point in-phase solutions {δ_s, α_s, g_01s, g_10s} = {18.8, 0.1, 0.073, 0.027}. We find that all the eigenvalues are real and {λ_δ, λ_α} < 1, indicating stability. However, {λ_g01, λ_g10} = 1. These eigenvalues suggest that while the in-phase solution is stable for static synapses, with STDP the system may be marginally stable or unstable. In order to determine if the system was marginally stable at the in-phase solutions, we examined the Jordan form of the linearized state matrix and observed that the Jordan blocks corresponding to unit eigenvalues were scalar. This served as theoretical confirmation for local stability of the system with STDP present.

Additionally, we numerically analyzed the sensitivity of the in-phase solutions reached through STDP. We tested the sensitivity of the solutions by varying the initial conditions of the state variables in Equation (15). In Figure 9A we show the evolution of the 1:1 in-phase map period of the 2-MCI network: 〈T_0∕1〉 = δ_S for different initial values of η₀ and H = 0. We see that for η₀ = 0 and even η₀ = 40, 〈T_0∕1〉 evolves to the stable in-phase period of 18.9 ms, which is the exact in-phase period predicted by the 2-MCI simulations. The map stably evolves to the in-phase solution even for an initial value that is 200% greater than the solution point. We also confirm that 〈T_0∕1〉 evolves to the in-phase period for different initial spike time lag (δ₀) values of {0, 4, 8, 12} ms. This is illustrated in Figure 9B. These numerical analyses confirm that the in-phase state of the 2-MCI system is indeed locally stable.

FIGURE 9

Figure 9. Numerical stability analysis of the 2-MCI map. (A) The evolution of the map's predicted synchrony period 〈T_0∕1〉 for different initial conditions of η₀. (B) The evolution of the map's synchrony period 〈T_0∕1〉 for different initial conditions of δ.

3.3. STDP Induced Emergent Synchronization in 100-MCI Network

In this section we examine how STDP allows larger MCI networks to synchronize. We consider a homogeneous all-to-all coupled network of 100 neurons with temporal heterogeneity H = 10. In the absence of STDP, Figure 10A shows a raster of the spike times (magenta) when the network fails to synchronize. On enabling STDP, the network is able to reach perfect in-phase synchronization as shown by the black rasters in Figure 10A. In Figure 10B we show the fraction of η values before (magenta) and after STDP (black). We can see that STDP evolves the η-values toward a slightly positive-skewed Gaussian distribution with a mean of -20. As we have seen in the previous sections, STDP tends to evolve the synapses such that η decreases. As η decreases, the spike times become closer to each other (within the 1:1 synchronization domain). When η is sufficiently negative, the spike time differences are almost zero. This causes the STDP to suddenly stop evolving, causing the distribution of η to appear positively skewed.

FIGURE 10

Figure 10. STDP based synchronization in a 100-MCI network. (A) Shows a raster of the spike times with STDP (black) and without STDP (magenta). (B) The corresponding distribution of η-values with STDP (black) and without STDP (magenta).

In Figure 11, we show how increasing H affects the ability of the 100 neuron MCI network to exhibit in-phase synchronization. In the absence of STDP, the synchronization measure S very quickly drops with increased H. However, when STDP is enabled, the network is able to sustain in-phase synchronization for H < 20. For H > 20, the network's synchronization measure gradually decreases.

FIGURE 11

Figure 11. 100-MCI Synchrony metric. Here we measure the measured synchrony of the 100-MCI network for a range of 0 ≤ H ≤ 50. In-phase synchronization is significantly enhanced for a wide range of 0 ≤ H ≤ 18 when STDP is enabled (magenta).

We also look at the other network metrics that we defined in the methods section. In Figure 12, we show the network link imbalance and neuronal strengths. The link imbalance measures the difference in synaptic strengths of coupled neurons in the network. It is expressed as a matrix L_ij = a_ijg_ij − a_jig_ji, where i and j refer to the pre-synaptic and post-synaptic neurons, respectively, and [a_ij] ∈ [0, 1] refers to the connectivity between neuron i and j. For a network with homogenous synaptic strengths, L is a zero matrix. When the network evolves under STDP to in-phase synchronization, we notice in Figure 12A that the network imbalance shifts from a zero matrix to a skew symmetric matrix whose main diagonal is zero (indicating no self-synapses). For i < j, L_ij > 0 and for i > j, L_ij < 0. This indicates that STDP evolves such that the strength of synapses from slower neurons onto faster neurons are larger in value. This observation is further reinforced by the neuronal strength metric ($G_i=\sum _j{a}_{ij}{g}_{ij}$) shown in Figure 12B. Here, the sum of the strengths of outgoing synapses for every neuron is presented. Without STDP, all neurons project equal synaptic influences on each other. With STDP, as consistent with our earlier observations, the slowest neuron number 0 has the largest synaptic influence on the rest of the network while the fastest firing neuron number 100 has the weakest synaptic influence on other neurons. This shows that with STDP, the slower neurons can control the firing rate of the faster neurons thus increasing the chances of attaining synchronization. Our 100-MCI network investigations suggest that STDP is an effective mechanism that allows in-phase synchronization for a wide range of neuron firing rate heterogeneity. It also suggests that STDP induced synchronization may not be significantly affected by the scale of the MCI network.

FIGURE 12

Figure 12. Network metrics. (A) Link imbalance of the 100-MCI network when H = 10. The synaptic strength differences between pairs of neurons are color coded. In general synapses originating from slower neurons (lower numbers on x-axis) have a positive network imbalance value. (B) The neuron strength metric measures the total outgoing synaptic strength for each neuron. With STDP, neuron strength linearly decreases in the direction of slowest firing neuron to the fastest neuron, indicating greater synaptic influences of slower neurons.

4. Discussion

In this study we investigated if STDP in mutually coupled interneuronal networks can induce stable in-phase synchronization at gamma frequencies. We first investigated the domains of synchronization for a 2-MCI network in terms of structural and temporal heterogeneity, given by {η, H}. We identified the 1:1, in-phase and higher order synchronization domains for a static 2-MCI network. We noted that in general, as H increased, η had to decrease in order to maintain 1:1 and in-phase synchronization. The decrease in η corresponds to the case of g₀₁ increasing while g₁₀ decreases. This in turn, indicated that, for 1:1 and in-phase synchrony, the slower neuron 0, suppresses the faster neuron 1, such that their periods match. In the case of higher order synchronization such as 2:1 synchronization, the faster neuron fires more than once for every single time the slower neuron fires. Here, the period of the slower neuron is significantly increased due to suppression by the faster neuron.

Previous work by White et al. (1998) has shown that synchronization in MCI networks are susceptible to even slight increase in firing rate heterogeneity. Our group had shown that STDP was one way by which stable synchronous oscillations could be achieved in UCI and 2-MCI networks (Talathi et al., 2008). This was based on the findings of Haas et al. (2006), that reported STDP was present in inhibitory synapses in the second layer of the entorhinal cortex. Experimental measurements provided the STDP learning rule (Haas et al., 2006) that was used in computational experiments by Talathi et al. (2008). We demonstrated that in the presence of firing rate heterogeneity, STDP resulted in significant enhancement of network synchronization. A STRC based theoretical map was derived to verify the stability of UCI in-phase synchronization.

In other works, STRC's have been used to calculate future spike times based on periodic synaptic perturbations from other neurons. Most of these studies assumed that the STRC's displayed weak second order components. For example in the work done by Talathi et al. (2008), due to weak second and higher order STRC components, the total spike time change for periodic synaptic perturbations was represented as Φ_i∞(δt, g) = Φ_i1(δt, g) + Φ_i2(δt, g). Where, δt was the time of perturbation after every spike in each firing cycle. In the course of our investigations we observed that a fast firing Wang-Buzsaki neuron (Wang and Buzsáki, 1996) with inhibitory synaptic input, produced STRC's with a significant second order component. This in turn indicated that calculating spike time shifts using the Φ_i∞ term would lead to incorrect predictions.

Recently, Talathi et al. (2009) demonstrated an empirical approach for predicting spike times when higher order STRC's were present. However, this method was specific for shunting synapses. The approach used in this paper is more generic. We considered cases where a 2-UCI network exhibited synchrony for a certain set of parameters. Then the error between predicted spike times and actual spike times were noted, and an approximately linear relationship between perturbation time δ_t and error in spike-times $\overline{{\delta }_t}$ was observed. Using this relationship we derived a correction factor for the perturbation time, which when incorporated into STRC calculation, allowed the term Φ_i∞ to more precisely describe spike time shifts.

In this work, we performed an in-depth analysis on the synchronization of 2-MCI network and eventually, a 100-MCI network. We utilized the same STDP learning rule form as used by Talathi et al. (2008). We also use a more biophysically relevant fast spiking interneuron model of Parvalbumin-expressing basket cells (Wang and Buzsáki, 1996). As in our previous work, the STDP mechanism implemented makes two assumptions: (1) Only the spike time difference between neighboring pairs of spikes was considered. (2) The effects of STDP modulation were linearly summed. These assumptions allowed us to obtain an analytical form of the STDP rule.

We began our study on the 2-MCI network by identifying the synchronization domain in terms of the structural and firing rate heterogeneity parameters {η, H}, respectively. We showed that decreasing η could compensate for increasing H (up to an extent), in terms of maintaining 1:1 and in-phase synchronization. We also identified higher order synchronization domains for the 2-MCI network. The synchronization manifold of the 2-MCI network provided us with an understanding of the direction in which η would have to evolve in a 2-MCI network, in order to attain 1:1 or in-phase synchronization, as H increased.

In the presence of STDP, for η₀ = 0, the 2-MCI network was found to evolve to in-phase synchronization for a significant range of H < 25. We determined that for H > 25, STDP had a higher probability to evolve the network toward higher order synchronization. Regardless, STDP enhanced the domain of synchronization not only for 1:1 synchronization, but also for higher order synchronization. This also indicated that depending on the initial conditions of η₀, STDP could variably evolve the network to different orders of synchrony.

In order to determine if the 1:1 and in-phase synchronization states of the network were stable, we mathematically analyzed the 2-MCI network and derived a STRC based map for predicting 1:1 and in-phase synchronization. We first validated the map by ensuring that the predicted spike time differences ({δ, α}) were the same as that of the static 2-MCI simulations, within the 1:1 synchronization domain. We then observed the map's evolution of η with STDP for an initial η₀ = 0 and H = 10. We found that the predicted value of η for in-phase synchronization, was around 44 which was close to the 2-MCI STDP evolution of η = 40. The slight discrepancy in the values of η most likely arose from the linear approximation of the STRC correction term $\overline{\delta t}$ as a function of δt_s. Having validated that the map behaves closely to the 2-MCI simulations, we performed linear stability analysis on the discrete mathematical map expressed in Equation (15). Our findings indicated that the eigenvalues for the spike times {δ, α} were stable (< 1). However, the eigenvalues corresponding to the STDP evolution of {g₀₁, g₁₀} were marginally stable. To prove that system was at least locally stable, we performed further numerical analyses. Specifically, we examined the evolution of the map with STDP for different initial conditions of δ and η. We observed that the map's predicted 1:1 synchronization period (〈T_0∕1〉 ≈ δ) and confirmed that it evolved to the period of the 2-MCI simulation in-phase synchronization, for a range of initial conditions. In fact η had to be increased significantly (η ≥ 20) beyond the stable value η = −40 in order for the map to fail to evolve to the in-phase synchronization regime.

We next examined how the scale of MCI networks affects STDP induced synchronization. We constructed and simulated a 100-MCI structurally homogeneous network (η₀ = 0) with and without STDP for H = 10. In accordance with our previous results, we found that without STDP, the structurally homogeneous network was completely desynchronized. On enabling STDP, the 100-MCI network was able to exhibit in-phase synchronization for a significantly increased range of temporal heterogeneity H < 18. It must be noted that as η is only defined for any given pair of mutually coupled interneurons, in a larger 100-MCI network (without self synapses) the total number of synaptic pairs are 4950, and each pair has an associated η that describes the structural heterogeneity. We observed that with STDP, the distribution of η evolved toward a normal distribution with a slight positive skew. This small positive tail of η-values indicates that the distribution of η evolved toward the negative direction and ceased evolving when in-phase synchronization was reached.

All of our results suggest that STDP is a very viable mechanism for the formation of interneuronal gamma (ING) oscillations. The original firing rates of the free running interneurons (H = 0) are approximately 60 Hz. In the 2-MCI network, we observed that the firing rates were in the ranges of 52 and 53 Hz, for 1:1 and in-phase synchronization, respectively. The 100-MCI network was also able to sustain 50 Hz in-phase synchronization. These frequencies are well within the gamma band and prove that an ING mechanism in conjunction with interneuronal STDP might be a strong candidate for stable gamma oscillations in cortical and even hippocampal regions (Bragin et al., 1995). A future research endeavor would be the inclusion of fast firing excitatory pyramidal cells network along with our MCI model to confirm that the ING mechanism can indeed drive pyramidal cells to synchronize at gamma frequencies. Additionally, it would be of great interest to see if and how STDP can induce 1:1 and in-phase synchronization in non-MCI networks with more random architectures. In conclusion, this work provides further evidence that for significant heterogeneity in firing rates, STDP may be a viable alternative to other synchronization mechanisms such as gap junctions (Traub et al., 2001).

Author Contributions

ST, PK conceived the initial design of the work with additional input from SR. SR performed all the simulation experiments and analysis. SR wrote the manuscript with significant inputs from PK and ST.

Funding

This research was funded by startup funds to ST from the Department of Pediatrics at the University of Florida; the intramural grant on Computational Biology at the University of Florida; and the Wilder Center of Excellence for Epilepsy Research and the Children's Miracle Network. PK was partially supported by the Eckis Professor Endowment at the University of Florida.

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

Abarbanel, H., Gibb, L., Huerta, R., and Rabinovich, M. (2003). Biophysical model of synaptic plasticity dynamics. Biol. Cybern. 89, 214–226. doi: 10.1007/s00422-003-0422-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Acker, C., Kopell, N., and White, J. (2004). Synchronization of strongly coupled excitatory neurons: relating network behavior to biophysics. J. Comp. Neurosci. 15, 71–90. doi: 10.1023/A:1024474819512

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartos, M., Vida, I., Frotscher, M., Meyer, A., Monyer, H., Geiger, J., et al. (2002). Fast synaptic inhibition promotes synchronized gamma oscillations in hippocampal interneuron networks. Proc. Natl. Acad. Sci. U.S.A. 99, 13222–13227. doi: 10.1073/pnas.192233099

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartos, M., Vida, I., and Jonas, P. (2007). Synaptic mechanisms of synchronized gamma oscillations in inhibitory interneuron networks. Nat. Rev. Neurosci. 8, 45–56. doi: 10.1038/nrn2044

PubMed Abstract | CrossRef Full Text | Google Scholar

Bi, G. Q., and Poo, M. M. (1998). Synaptic modifications in cultured hippocampal neurons: dependence on spike timing, synaptic strength, and postsynaptic cell type. J. Neurosci. 18, 10464–10472.

PubMed Abstract | Google Scholar

Bi, G. Q., and Wang, H. X. (2002). Temporal asymmetry in spike timing-dependent synaptic plasticity. Physiol. Behav. 77, 551–555. doi: 10.1016/S0031-9384(02)00933-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Bragin, A., Jandó, G., Nádasdy, Z., Hetke, J., Wise, K., and Buzsáki, G. (1995). Gamma (40–100 hz) oscillation in the hippocampus of the behaving rat. J. Neurosci. 15(1 Pt 1), 47–60.

PubMed Abstract | Google Scholar

Caporale, N., and Dan, Y. (2008). Spike timing-dependent plasticity: a hebbian learning rule. Annu. Rev. Neurosci. 31, 25–46. doi: 10.1146/annurev.neuro.31.060407.125639

PubMed Abstract | CrossRef Full Text | Google Scholar

Ermentrout, B. (1996). Type 1 membranes, phase resetting curves and synchrony. Neural Comput. 8, 979–1001.

PubMed Abstract | Google Scholar

Froemke, R. C., and Dan, Y. (2002). Spike-timing-dependent synaptic modification induced by natural spike trains. Nature 416, 433–438. doi: 10.1038/416433a

PubMed Abstract | CrossRef Full Text | Google Scholar

Froemke, R. C., Tsay, I. A., Raad, M., Long, J. D., and Dan, Y. (2006). Contribution of individual spikes in burst-induced long-term synaptic modification. J. Neurophysiol. 95, 1620–1629. doi: 10.1152/jn.00910.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Golomb, D., and Rinzel, J. (1993). Dynamics of globally coupled inhibitory neurons with heterogeneity. Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 48, 4810–4814.

PubMed Abstract | Google Scholar

Haas, J. S., Nowotny, T., and Abarbanel, H. D. I. (2006). Spike-timing-dependent plasticity of inhibitory synapses in the entorhinal cortex. J. Neurophysiol. 96, 3305–3313. doi: 10.1152/jn.00551.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Hansel, D., and Sompolinsky, H. (1992). Synchronization and computation in a chaotic neural network. Phys. Rev. Lett. 68, 718–721.

PubMed Abstract | Google Scholar

Kopell, N., and Ermentrout, B. (2004). Chemical and electrical synapses perform complementary roles in the synchronization of interneuronal networks. Proc. Natl. Acad. Sci. U.S.A. 101, 15482–15487. doi: 10.1073/pnas.0406343101

PubMed Abstract | CrossRef Full Text | Google Scholar

Kudela, P., Franaszczuk, P. J., and Bergey, G. K. (2003). Changing excitation and inhibition in simulated neural networks: effects on induced bursting behavior. Biol. Cybern. 88, 276–285. doi: 10.1007/s00422-002-0381-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Lytton, W. W., and Sejnowski, T. J. (1991). Simulations of cortical pyramidal neurons synchronized by inhibitory interneurons. J. Neurophysiol. 66, 1059–1079.

PubMed Abstract | Google Scholar

Miltner, W. H., Braun, C., Arnold, M., Witte, H., and Taub, E. (1999). Coherence of gamma-band eeg activity as a basis for associative learning. Nature 397, 434–436.

PubMed Abstract | Google Scholar

Oprisan, S., and Canavier, C. (2001). Stability analysis of ring of pulse coupled oscillators: the effect of phase resetting in the second cycle after the pulse is important at synchrony and for long pulses. Differ. Equat. Dynam. Syst. 9, 243–258.

Oprisan, S., Prinz, A., and Canavier, C. (2004). Phase resetting and phase locking in hybrid circuits of one model and one biological neuron. Biophys. J. 87, 2283–2298. doi: 10.1529/biophysj.104.046193

PubMed Abstract | CrossRef Full Text | Google Scholar

Pulvermller, F., Lutzenberger, W., Preissl, H., and Birbaumer, N. (1995). Spectral responses in the gamma-band: physiological signs of higher cognitive processes? Neuroreport 6, 2059–2064.

PubMed Abstract | Google Scholar

Ritz, R., and Sejnowski, T. J. (1997). Synchronous oscillatory activity in sensory systems: new vistas on mechanisms. Curr. Opin. Neurobiol. 7, 536–546.

PubMed Abstract | Google Scholar

Singer, W., and Gray, C. M. (1995). Visual feature integration and the temporal correlation hypothesis. Annu. Rev. Neurosci. 18, 555–586.

PubMed Abstract | Google Scholar

Talathi, S., Hwang, D., and Ditto, W. (2008). Spike timing dependent plasticity promotes synchrony of inhibitory networks in the presence of heterogeneity. J. Comp. Neurosci. 25, 262–281. doi: 10.1007/s10827-008-0077-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Talathi, S., Hwang, D., Miliotis, A., Carney, P., and Ditto, W. (2009). Predicting synchrony in heterogeneous pulse coupled oscillators. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 80:021908. doi: 10.1103/physreve.80.021908

PubMed Abstract | CrossRef Full Text | Google Scholar

Talathi, S. S., Hwang, D.-U., Carney, P. R., and Ditto, W. L. (2010). Synchrony with shunting inhibition in a feedforward inhibitory network. J. Comput. Neurosci. 28, 305–321. doi: 10.1007/s10827-009-0210-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Traub, R. D., Spruston, N., Soltesz, I., Konnerth, A., Whittington, M. A., and Jefferys, G. R. (1998). Gamma-frequency oscillations: a neuronal population phenomenon, regulated by synaptic and intrinsic cellular processes, and inducing synaptic plasticity. Prog. Neurobiol. 55, 563–575.

PubMed Abstract | Google Scholar

Traub, R. D., Kopell, N., Bibbig, A., Buhl, E. H., LeBeau, F. E., and Whittington, M. A. (2001). Gap junctions between interneuron dendrites can enhance synchrony of gamma oscillations in distributed networks. J. Neurosci. 21, 9478–9486.

PubMed Abstract | Google Scholar

van Vreeswijk, C., Abbott, L., and Ermentrout, B. (1994). When inhibition and not excitation synchronizes neural firing. J. Comp. Neurosci. 1, 313–321.

PubMed Abstract | Google Scholar

Wang, X., and Buzsáki, G. (1996). Gamma oscillation by synaptic inhibition in a hippocampal interneuronal network model. J. Neurosci. 16, 6402–6413.

PubMed Abstract | Google Scholar

Wang, X., and Rinzel, J. (1992). Alternating and synchronous rhythms in reciprocally inhibitory model neurons. Neural Comput. 4, 84–97.

Google Scholar

White, A., Chow, C., Ritt, J., Trevino, C., and Kopell, N. (1998). Synchronization and oscillatory dynamics in heterogeneous, mutually inhibited neurons. J. Comp. Neurosci. 5, 5–16.

PubMed Abstract | Google Scholar

Whittington, M. A., Traub, R. D., Faulkner, H. J., Stanford, I. M., and Jefferys, J. G. (1997). Recurrent excitatory postsynaptic potentials induced by synchronized fast cortical oscillations. Proc. Natl. Acad. Sci. U.S.A. 94, 12198–12203.

PubMed Abstract | Google Scholar