## Effect of Phase Response Curve Shape and Synaptic Driving Force on Synchronization of Coupled Neuronal Oscillators.

Monday, June 19, 2017

Ramana Dodla, Charles J. Wilson Neural Comput. 29(7):1769-1814.

The role of the phase response curve (PRC) shape on the synchrony of synaptically coupled oscillating neurons is examined. If the PRC is independent of the phase, because of the synaptic form of the coupling, synchrony is found to be stable for both excitatory and inhibitory coupling at all rates whereas the antisynchrony becomes stable at low rates. A faster synaptic rise helps extend the stability of antisynchrony to higher rates. If the PRC is not constant, but has a profile like that of a leaky integrates-and-fire model, then, in contrast to the earlier reports that did not include the voltage effects, mutual excitation could lead to stable synchrony provided the synaptic reversal potential is below the voltage level the neuron would have reached in the absence of the interaction and threshold reset. This level is controlled by the applied current and the leakage parameters. Such synchrony is contingent on significant phase response (that would result, for example, by a sharp PRC jump) occurring during the synaptic rising phase. The rising phase however does not contribute significantly if it occurs before the voltage spike reaches its peak. Then a stable near-synchronous state can still exist between type-1 PRC neurons if the PRC shows a left-skewness in its shape. These results are examined comprehensively using perfect integrate-and-fire, leaky integrate-and-fire and skewed PRC shapes under the assumption of the weakly coupled oscillator theory applied to synaptically coupled neuron models.

Figure 1: Numerical simulation of a network of two “perfect” integrate-and-fire neuron models. Coupling is excitatory, but results are qualitatively similar when the coupling is inhibitory. (a) The two voltage time courses at I0 = 0.1 approaching a synchronous state. (b) The same initial conditions as in (a) leading to antisynchronous state, but a more closely spaced initial conditions leading to a synchronous state when I0 is lowered. Synaptic conductance is an alpha function with α = 1/3. Vreset = 0, Vth = 1. Esyn = 2.

Figure 2: Analytical prediction of synchrony and anti-synchrony between two perfect integrate-and-fire neuron models. (a) The interaction function H(φ) and the alpha-function synaptic conductance at an Esyn that makes the coupling excitatory. The PRC of the model is simply 1/Iapp. (b) The growth function G(φ) at three different Iapp. G(φ), the phase-locked states and their stability are independent of Esyn, and thus the panels (b-d) remain unaltered between excitatory and inhibitory couplings. (c) The steady phase differences as a function of the frequency displaying stable synchrony (filled circles resembling thick lines at ordinate values of 0 and T) at all frequencies, stable antisynchrony (filled circles at T/2) at low frequencies. The unfilled circles represent unstable steady state phase differences. (d) Parameter display of synchrony (shaded, filling all the space) and antisynchrony (hatched) in the space of synaptic time constant τ (= 1/α) and firing rate (= 1/T ). Unit of ms is assigned to time for convenience.

Figure 3: Analytical prediction of antisynchrony between two perfect integrate-and-fire neuron models using double exponential synaptic function in the plane of firing rate and synaptic decay time. The region to the left of each curve (the hatched region) marks the stability of antisynchrony for the corresponding synaptic rise time. Synchrony is stable in the entire parameter space.

Figure 4: Simulation of a mutually excitatory network of two leaky integrate and fire neurons whose voltages are modeled according to those of the Wang-Buzsáki model. (a) The voltage time courses [V1(t), V2(t)] evolving in time and approaching a synchronous oscillatory state at a firing frequency of 100 Hz. Spikes are not part of the model, and hence are not drawn. On the right, the effect of the coupling on V1(t) is shown - it is excitatory. (b) Time courses as in (a) but at 50 Hz. Two sets of initial conditions lead to two different states: antisynchrony for widely separated initial conditions, and synchrony for closely separated initial conditions. Right: The difference of successive spike times of the neurons normalized to the instantaneous period of one of the neurons (i.e the phase difference) is shown for the two sets of initial conditions. (c) As in (b) but at 10 Hz. Synchrony and anti synchrony are both unstable, but a near synchronous state that is very close to the synchronous state is stable. For all the simulations The parameters are GL = 0.01 mS/cm2, EL = 0 mV, Cm = 1 μF/cm2, Vreset = −100 mV, Vth = −49.5635 mV. Esyn = 10 mV, α = 1/3 ms−1, and g = 0.04 mS/cm2.

Figure 5: Simulation of a pair of mutually excitatory LIF model neurons with the conventional parameters without units: Cm = 1, GL = 1, EL = 0, Vreset = 0, and Vth = 1. (a) Voltage time courses of the two neuron models exhibiting antisynchrony, non-zero phase-locked state, and synchrony at three different levels of the applied current. (b) The phase difference of the consecutive spike times of the illustrations in (a) displaying the evolution toward steady states. The other parameters are Esyn = 2 that is above the threshold, and α = 6 that determines the shape of the alpha-function synaptic conductance. The membrane effect is included in the computation.

Figure 6: Analytical prediction of the dynamics of a network of two coupled leaky integrate-and-fire neuron models whose voltages are shaped according to those of the Wang-Buzáski model. (a) The voltage time course of the each uncoupled LIF model (thick lines) compared to that of Wang-Buzs ́aki model (thin lines). (b) The PRC of the LIF model. (c) The firing rate as a function of the applied current I0. (d,e)Mutual excitation. (d) Growth functions at three levels of I0. I0 = 4.3 μA/cm2 for 100 Hz, 1.7825 μA/cm2 for 50 Hz, and 0.53 μA/cm2 for 25 Hz. (e) Stable (filled circles) and unstable (unfilled circles) steady state phase differences with frequency. A range of firing rates (or I0) can result in stable synchrony. (f,g) Mutual inhibition. (f) Growth functions are illustrated at three frequencies. (g) Steady state stability of synchrony and other states with I0.a

Figure 7: Analytical prediction of the dynamics of a pair of mutually excitatory LIF model neurons that use the conventional parameters as those set in Fig. 5. (a) Voltage time course of each uncoupled model neuron, (b) phase response curve of each neuron, (c) firing rate of each uncoupled neuron with steady input current. If time were measured in ms, then the firing rate would be measured in kHz, and achieving firing rates of the order of 10’s of Hz would become an extremely sensitive function of I0. (d) Growth functions at Esyn = 2, α = 1/3 displaying stable synchrony as I0 is increased across (Esyn − EL)GL. (e) Steady state phase differences with I0. Filled circles indicate stability, and unfilled instability. (f) Steady state phase differences with α at I0 = 1.5, and (g) same as in (f) but with I0 = 2.5. The results in (f) are similar to those of Fig. 1 of Van Vreeswijk et al. (Van Vreeswijk et al., 1994).

Figure 8: Regions of stability of synchrony and antisynchrony for the LIF model with the conventional parameters. (a) Stability region in I0 vs. α plane as the reversal potential is decreased from begin mutual excitatory to being mutual inhibitory. The shaded region is the region of stable synchrony, hatched region is the region of stable antisynchrony. Multiple regions of bistability, as well as non-zero phase-locked states (white space) also exist. (b) Stability regions in I0 and Esyn space at different levels of α. (c) Stability regions in α and Esyn space at different levels of the applied current. The model parameters are as those in Fig. 7.

Figure 9: Near-synchrony in Wang-Buzs ́aki excitatory network. (a-c) Failure of synchrony despite sharp rise of the PRC near small phases, or delay in conductance. ε = 5, I0 = 0.211 μA/cm2. (a) V (t), the synaptic conductance of the Wang-Buzs ́aki model illustrated with and without delay, and PRC profile at small phases. (b) The growth function of two mutually excitatory neurons without and with (d = 1 ms) delay. (c) Spike times of a network of 100 coupled neurons without synaptic delay failing to exhibit synchrony. g ̄ = 0.02 mS/cm2, Esyn = 0 mV. The initial conditions of are chosen randomly on the limit cycle orbit of the uncoupled model. (d-f) Spike times of the Wang-Buzs ́aki network as the skewness of the PRC is altered while keeping the firing rate at 10 Hz (by adjusting Iapp) and d = 0. Corresponding to these three simulations, the PRCs are shown in (g), and the growth functions for two coupled oscillators are shown in (h). (i) The stable (filled circle) and unstable (unfilled circle) phase-locked solutions as a function of ε. d = 0.

Figure 10: Effect of PRC skewness on synchrony. (a) PRC shapes given by Eq. 39 as n is increased. (b) The growth function for two mutually excitatory neurons getting affected by the period and in particular the antisynchrony becoming unstable as T is increased even when the PRC has a small skewness. (c) Change of stability of synchrony, antisynchrony, and the non-zero phase-locked states with the rate for excitatory network when n = 1. (d) Same as in (c) for inhibitory network. (e,f) Similar to (c) and (d) but for large skewness (n = 10), depicting only the movement at low rates. In (e) the non-zero phase-locked state is in close proximity to the unstable synchronous state (which is weakly repelling) at low rates making it a near-synchronous state. (g) The growth functions for excitatory network revealing the non-zero phase locked state surging toward the synchronous state as the rate is lowered at large skewness (n = 10). (h) Movement of steady states of the excitatory network as the skewness is increased at 10 Hz. (i) Steady states displayed in n vs. T space for the excitatory network.

Figure 11: Subdued PRC response near early phases causes near-synchrony. (a, top) PRCs as in Eq. 39 for n = 0, 0.5, 1 but are reset to the case of n = 0 for phases above T /2. Period of oscillation is 100 ms. The peak levels of the PRCs are normalized to unity. (a, bottom) Corresponding interaction functions for mutual excitation using double exponential synaptic function with τd = 3 ms and τr = 0.1 ms. The interaction functions are computed without including the voltage effect. (c) The detail of the H(φ) and H(−φ) near the early phases. (d) Growth functions corresponding to the PRCs in (a). (d-f) As in (a-c) but resetting the PRCs to the case of n = 0 for phases below T/2.