Phase resetting curves allow for simple and accurate prediction of robust N:1 phase locking for strongly coupled neural oscillators

Abstract

Existence and stability criteria for harmonic locking modes were derived for two reciprocally pulse coupled oscillators based on their first and second order phase resetting curves. Our theoretical methods are general in the sense that no assumptions about the strength of coupling, type of synaptic coupling, and model are made. These methods were then tested using two reciprocally inhibitory Wang and Buzsáki model neurons. The existence of bands of 2:1, 3:1, 4:1, and 5:1 phase locking in the relative frequency parameter space was predicted correctly, as was the phase of the slow neuron's spike within the cycle of the fast neuron in which it occurred. For weak coupling the bands are very narrow, but strong coupling broadens the bands. The predictions of the pulse coupled method agreed with weak coupling methods in the weak coupling regime, but extended predictability into the strong coupling regime. We show that our prediction method generalizes to pairs of neural oscillators coupled through excitatory synapses, and to networks of multiple oscillatory neurons. The main limitation of the method is the central assumption that the effect of each input dies out before the next input is received.

Figure 1

(A) PRC generation. The top trace shows the membrane potential during PRC generation. The horizontal dotted line indicates zero mV. The parameters were as in the Methods section except that I_app = 0.5 μA/cm² for both neurons and g_syn = 0.15 mS/cm². The lower trace shows the perturbation in the dimensionless activation variable (s) for the synaptic conductance applied at a phase of ϕ = 0.5. The intrinsic period is P₀, the cycle containing the perturbation is P₁, and subsequent cycles are numbered sequentially. The stimulus time, ts, is the time between the previous spike and stimulus onset. (B) Typical phase resetting curve. The first order resetting is the solid line and the second order resetting is the dashed line. Third order resetting was not visible on this scale and is therefore negligible. The same parameters as in A were used.

Figure 2

Illustration of averaged weak coupling. I_app,F = 1.842 μA/cm² for the fast (F) neuron (solid line) and I_app,S = 0.77 μA/cm² for the slow (S) neuron (dashed line). These parameter values are for two oscillators where neuron S has a free running period which is double the period of neuron F. (A) The free-running membrane potential waveforms for the fast (solid line) neuron and the slow (dashed line) neuron. This potential was used in two ways: first to obtain the postsynaptic voltage V(ϕ_j) to calculate the postsynaptic current i(ϕ_i, ϕ_j) = g_syns(ϕ_i) {V(ϕ_j) − E_syn} and second to calculate the presynaptic voltage V(ϕ_i) to drive the presynaptic activation s(ϕ_i). (B) The presynaptic activation calculated as described in A. Notice that the solid black trace in A drives the dotted trace in B and vice versa. (C) The adjoint (iPRC) is computed for each neuron. (D) The coupling functions computed using the averaging method for weakly coupled oscillators for each neuron are shown here. This is the only step in which the value of g_syn is taken into account. The perturbations in the synaptic current are convolved with the iPRC over the network period to calculate these functions. They were numerically calculated using XPPAUT.

Figure 3

Mode prediction for weak averaged coupling. The solid black curve is the effective coupling function K(ϕ) for the system composed of the two neurons from Fig. 3 whose intrinsic periods have a ratio of exactly 2:1 for weak coupling with g_syn = 0.01 mS/cm². The lack of zero crossing shows that they cannot lock unless the intrinsic frequency of one neuron is modified. The dashed curves are the coupling function Q(ϕ) computed with additional current to the fast neuron to change its intrinsic frequency sufficiently to enable a locking. The range of current added to the fast neuron (−0.026 μA/cm² for the top dashed curve and −0.0385 μA/cm² for the lower dashed curve) that enables these curves to have a zero crossing of the x axis determines the lower and upper limits of the intrinsic frequency of the fast neuron that allow a locking. The region inside the dashed curves is the existence region for this weak coupling strength. The upward zero crossings indicate the phases (separated by ∼0.5) at which the slow neuron receives an input from the fast neuron in a stable 2:1 locking for the coupled system.

Figure 4

Mode prediction for pulsatile coupling. (A) Assumed firing pattern for N:1 locking. The faster neuron (F) fires N times for each time that the slower neuron (S) fires. Firing times are indicated by thick vertical bars. The input phases ϕ are defined on the upper and lower traces for dotted lines corresponding to a vertical bar on the partner trace indicating the time of a spike in the partner. Only the first and last spikes of the fast neuron within each cycle of the slow neuron are shown, the remaining spikes are indicated by the dots. The stimulus (ts) and recovery (tr) intervals are the time elapsed between vertical dotted lines. Both phases and intervals are subscripted by the neuron and indexed by the cycle. The second subscript on the phases of the slow neuron indicates the jth input during the current cycle of the slow neuron (cycles m and m+ 1 are shown). The intervals highlighted in gray are defined in the text and used to construct a discrete map from ϕ_SN[m] to ϕ_SN[m + 1], the phases highlighted in gray. (B–D) show an example prediction for 2:1 modes at E_syn = 0.75 mV, g_syn = 0.25 mS/cm², I_app = 1.0 μA/cm², and ɛ = 0.241 μA/cm². (B) PRCs for the faster neuron (B1) and the slower neuron (B2). The phases for four modes identified in C are shown. The open symbols correspond to unstable modes in C and the asterisks to stable modes. (C) Graphical method for existence. The error is the difference between assumed and predicted values of ϕ_SN and is plotted as a function of the assumed ϕ_SN. The open circles indicate unstable fixed points and the asterisk indicates a stable fixed point. (D) Testing the predictions. The full set of differential equations for two coupled neurons was integrated and produced a 2:1 mode with the membrane potential for the fast neuron (dotted curves) and slow neuron (solid curves) exhibiting the predicted intervals (labeled horizontal bars). Technical note: The open triangle in C is not exactly a zero crossing, but rather was detected by checking for a change in the sign of the error between the leftmost point on the error curve and the rightmost point on the error curve. The leftmost point corresponds to a ϕ_S1 > 0 and the rightmost to a ϕ_SN < 1, so any near synchronous solution must lie between these two points and the algorithm declares synchronous solution to be at the endpoint with the least error.

Figure 5

Qualitative comparison of predicted and observed modes for pulsatile coupling method as heterogeneity and conductance strength are varied. (A) The predictions were generated using the methodology described in Fig. 4. The circles indicate stable predicted modes at each locking ratio. (B) The observations were generated by numerically solving the full system of differential equations that describe the two coupled Wang and Buzsáki neurons. I_app was held constant at 1 μA/cm² as ɛ and g_syn were varied, and other parameters were as in the Methods.

Figure 6

Different phasic relationships observed at 2:1 frequency ratios. (A) Near synchrony; 2:1 locking, near synchrony between the slow (solid line) neuron spike and one of the fast (dashed line) neuron spikes, slow neuron leads. (B) Near antiphase; 2:1 locking, slow neuron fires near the midpoint of two fast spikes in near antiphase. (C) Leapfrog. Complex locking (4:2) in which the firing order of the fast and slow neurons switches each time the slow neuron fires. (D) Near synchrony; 2:1 locking, near synchrony between the slow neuron spike and one of the fast neuron spikes, fast neuron leads. These points were sampled from the points in the previous figure at g_syn = 0.15 mS/cm².

Figure 7

Quantitative comparison of predicted and observed intervals for pulsatile coupling method. (A) The values for the successive intervals ts_F, tr_F1, and tr_F2 as defined in Fig. 2 are shown for selected parameter values from Fig. 5 (g_syn = 0.15 mS/cm², I_app = 1 μA/ cm²). Open circles indicate observations and crosses indicate predictions. (B) The smallest stimulus interval at g_syn = 0.2 mS/cm² from Fig. 5 is plotted against epsilon. Circles indicate predictions and crosses and the X indicate observed. Colors indicate locking mode as in Fig. 5. Black circles indicate incorrect predictions as described in the text.

Figure 8

Qualitative comparison of predicted and observed modes for pulsatile coupling method as both intrinsic frequencies are varied. (A) The predictions were generated using the methodology described in Fig. 4. The circles indicate stable predicted modes at each locking ratio. (B) The observations were generated by numerically solving the full system of differential equations that describe the two coupled Wang and Buzsáki neurons. The colors indicate the value of N. The solid lines correspond to the very narrow bands predicted and observed at g_syn = 0.05 mS/cm² whereas the bands composed of small open circles correspond to the broader bands predicted at g_syn = 0.15 mS/cm². The stronger inhibitory coupling enabled the observation of 4:1 and 5:1 locking, which were not observed at the weaker coupling strength in the range of intrinsic frequencies shown.

Figure 9

Qualitative comparison of predicted and observed modes for stronger pulsatile coupling as both intrinsic frequencies are varied. (A) The predictions were generated using the methodology described in Fig. 4. The circles indicate stable predicted modes at each locking ratio. (B) The observations were generated by numerically solving the full system of differential equations that describe the two coupled Wang and Buzsáki neurons. The colors indicate the value of N whereas the bands composed of small open circles correspond to the bands predicted at g_syn = 0.25 mS/cm².

Figure 10

Weak coupling predictions compared to observations from our method. Shown are the effective coupling function (solid black curve) for the coupled system with g_syn = 0.25 mS/cm². The system does not have a stable locking mode at the values of the applied current that give rise to a 2:1 mode in the uncoupled system. Changing the applied current value for the fast neuron shifts the coupling function down. For the range of values where there was a predicted and observed 2:1 locking in the system with our method (shaded region enclosed by the dashed curves), weak coupling method does not predict such a mode. This is due to a failure to satisfy the conditions required for the weak coupling method.

Figure 11

Accurate prediction of 2:1 locking band in two neuron network with excitation. The parameters were the same as in Figs. 8 and 9 except E_syn = 0 mV and g_syn = 0.04 mS/cm². Two to one lockings are indicated in red. For completeness, 1:1 lockings are shown in gray. (A) Predicted. Modes were predicted directly from the PRCs using the existence and stability criteria. (B) Observed. The solution to the differential equations were run until the transients dissipated (1 s) and locked modes were identified by the periodic repetition of intervals. The agreement between predicted and observed is excellent.

Figure 12

Prediction of population locking. (A) A four neuron network with two clusters of neurons with different intrinsic frequencies. Each cluster is indicated inside a box. The within cluster connections are inhibitory (solid circles) and the between cluster connections are excitatory (open triangles). Each cluster is collapsed into a single oscillator representing the cluster. (B) The example chosen shows 2:1 locking as predicted. The traces of the neurons in the fast (solid lines) and slow (dashed lines) cluster are indistinguishable.