Stabilizing synchrony by inhomogeneity

We show that for two weakly coupled identical neuronal oscillators with strictly positive phase resetting curve, isochronous synchrony can only be seen in the absence of noise and an arbitrarily weak noise can destroy entrainment and generate intermittent phase slips. Small inhomogeneity–mismatch in the intrinsic firing rate of the neurons–can stabilize the phase locking and lead to more precise relative spike timing of the two neurons. The results can explain how for a class of neuronal models, including leaky integrate-fire model, inhomogeneity can increase correlation of spike trains when the neurons are synaptically connected.

X 1 = F (X 1 ) + I 1 + g 12 G 12 (X 1 , X 2 ) + σξ 1 (t) X 2 = F (X 2 ) + I 2 + g 21 G 21 (X 2 , X 1 ) + σξ 2 (t) (S1) where X i is a N-dimensional state vector containing the membrane potential and gating variables. For example in the Hodgkin-Huxley (HH) model [S1], X = [V, m, h, n] T and in the Leaky Integrate-and-Fire (LIF) model [S2], X = V . F (X) defines the internal dynamics of the neuron i. G ij (X i , X j ) determines functional form of synaptic connection from neuron j to neuron i. For example in the case of pulse coupled synapses, as we used in this letter it would be G ij X i (t) , X j (t) = n δ t − t n j (S2) where t n j is the instant of n th spike of neuron j and δ (t) is the Dirac's delta function. g ij shows the synaptic weight from neuron j to neuron i which is scaled by a small factor ε so that the weakly coupled oscillators approximation is valid.
Each neuron receives suprathreshold constant current I i with mismatch ∆I = I 1 − I 2 as well as an independent Gaussian white noise characterized by its mean and auto-correlation, and its cross-correlation with the other neuron's input We assume that at the absence of synaptic connections and noise, dynamics of each isolated neuron with the input current I i given byẊ has a T-periodic limit cycle solution X 0 (t). Every point on the limit cycle X 0 can be characterized by a single phase variable θ. On the limit cycle X = X 0 (θ), which its inverse maps each point on the limit cycle to a unique phase, θ = X −1 0 (X) = Φ(X). The phase of each oscillator along the limit cycle is then where the relative phase φ i is a constant which is determined by the initial conditions [S3]. From the definition of phase (Eq. S6), the rate of the phase change on the limit cycle is So by means of the chain role we obtain the following useful relation: Where ∇ X Φ is the gradient of Φ with respect to the vector state variable X.
Using the notion of asymptotic phase and isochrons (See Refs. [S4] and [S3]) the phase of oscillation θ = Φ(X) can also be defined for all the points X in the neighborhood of the limit cycle. By this concept we can introduce Z(θ) = ∇ X Φ as the infinitesimal phase response curve (iPRC): If the neuron on the limit cycle which is in the state X(t), receives a small vector perturbation εU in a short time interval ∆t, immediately after the perturbation the state of the neuron jumps to X(t) + εU . So the resulted phase change is ∆θ ∇ X Φ (X(θ)) .εU (S9) Therefore in the presence of perturbation Eq. (S8) changes as We consider the coupling and noise terms in Eq. S1 as small perturbations on the limit cycleẊ i =F (X i ). We also assume that inherent frequency of neurons have a small difference of the order O(ε). (1)) we can reduce high dimensional Eqs. S1 to two scalar equations which determine the evolution of the phase [S5-S8] where Z (θ) is the infinitesimal phase response curve [S9] and ε∆ω = ω 2 − ω 1 . Assuming X i (t) X 0 (θ i (t)) and by the change of variable θ i (t) = ωt + φ i (t), the equations for the evolution of the relative phase of the oscillators read: We exploit the fact that ε is small to further reduce Eqs. S14. With a system of the forṁ Averaging theory [S8, S10, S11] states that in Eq. S15, x (t) can be replaced by its average over a periodx anḋ By applying averaging method on the Eqs. S14 we havė originates from averaging the noisy phase equations [S12], and without loss of generality we assumed that the phase is normalized so that 0 ≤ θ < T , i.e., ω = 1. By defining φ = φ 2 − φ 1 , we derive the following equation for the phase differencė where √ 2η (t) = ξ 2 − ξ 1 and η(t) is a Gaussian white noise with zero mean and unit variance.
, for QIF oscillator (top-left) and LIF oscillator (top-right) for different values of difference of the synaptic strengths. For the LIF oscillators for which PRC is not an even function, the coupling term can be non-zero for symmetric connections whereas for QIF oscillators with an even PRC the effective coupling is determined by the difference of two reciprocal synaptic constants. In the lower panel we have shown how increasing the mismatch in the firing rates of the neurons stabilizes the fixed point and then expands its basin of attraction.
We have used three model neurons in our study: (i) Canonical type-I oscillators with Z(φ) = 1 − cos(φ), (ii) LIF oscillators which is described byv = I − v for t ≤ v T h (= 1) and lim τ →0 + v(t + τ ) = 0 with where ω = 2π/[logI − log(I − 1)] [S13], and (iii) the conductance based Wang-Buzsaki (WB) model as described at the end of the supplementary material [S14]. First we focus on the deterministic case of Eq. S24 with D = 0. For QIF oscillator, Z (φ) = 1 − cos (φ) is an even function of φ. Therefore it reduces toφ where ∆g = g 21 −g 12 is the effective coupling constant. In this case the most effective coupling is that which maximizes ∆g. i.e., a unidirectional one and the symmetric connection leads to zero effective coupling ∆g = 0. But for LIF neurons with uneven PRC Eq. S24 takes the forṁ Note that for the oscillators with an oblique PRC, e.g. the LIF oscillators, the effective coupling term can be non-zero for symmetric connections (see Fig. S1). The fix points of general equation S24 with D = 0 are the intersections of horizontal line T ∆ω and the curve described by F (φ) = g 12 Z(φ) − g 21 Z(−φ).
For type-I phase oscillators we rewrite Langevin Eq. S24 as with Γ (φ) = ∆ω ∆g − 1 T + 1 T cos φ . Corresponding Fokker-Planck equation [S15-S17] for the phase difference distribution ρ(φ, t) is The stationary phase difference distribution satisfies with the solution where N is a normalization factor so that T 0 ρ (φ) dφ = 1, and α = Dσ 2 φ ∆g is the ratio of noise intensity to the coupling strength. The constant A can be determined by the periodicity condition of ρ 0 , that is, ρ 0 (0) = ρ 0 (T ). Therefore the final form of stationary solution is In Figure 2A we have shown the result of analytic solution for steady-stat phase difference distributions S33 and that of direct numerical integration of of phase differential equations S12. For solving Eq. S33, we have used double int function of MATLAB. In simulation, we integrate Eqs. S12 with Euler method and save spike times of each neuron. Then we have used hist function in MATLAB to plot ρ.
It has been shown that [S18] in the weak coupling and weak noise limit, the cross-correlogram (CC) and the phase difference probability distribution function, ρ(φ), are related by The most probable phase difference of spiking of the neurons (location of the peak of the PDF in Fig. 2A) can be determined by differentiation of ρ with respect to φ. Derivative of ρ with respect to φ is For and Eq. S35 reduces to and By using equations S39 and S33 we have We have used MATLAB function "int" to plot φ * versus ∆ω in Figure 2C using Eqs. S33, S41, and C max ; and ρ max versus ∆ω in Figure 2B by Eq. S40.

NEURONAL MODELS
In Figure 3 (left panel) we have plotted C(τ ) for two leaky integrate-fire (LIF) neurons described by . We integrated this equation for two pulse coupled neurons and calculated spike count for a time window of teh width 0.5mS. Parameters are selected in accordance with experimental data [S19] as τ m = 20 ms, v res = −74 mV , V th = −54 mV , v r = −60 mV , I = 25.0 + ∆I mV , g ij = 1 mV and D = 1.0. In figure 3 (right panel) we have used the conductance-based Wang-Buzsaki (WB) model [S14] which is widely used for simulation of Type-I neurons. The current balance equation for this model is where C is the membrane capacitance in µF/cm 2 , vi is the membrane voltage in mV . g N a , g K and g l are the maximum conductance per surface for the sodium, potassium and leak currents, respectively, and V N a , V K and V l are constants which determine the corresponding reversal potentials [S20]. m, and h are the activation variable for the inactivation variable for sodium current, respectively; and n is the activation variable for potassium current. The time course of the gating variables obey