Using phase resetting to predict 1:1 and 2:2 locking in two neuron networks in which firing order is not always preserved

Abstract

Our goal is to understand how nearly synchronous modes arise in heterogenous networks of neurons. In heterogenous networks, instead of exact synchrony, nearly synchronous modes arise, which include both 1:1 and 2:2 phase-locked modes. Existence and stability criteria for 2:2 phase-locked modes in reciprocally coupled two neuron circuits were derived based on the open loop phase resetting curve (PRC) without the assumption of weak coupling. The PRC for each component neuron was generated using the change in synaptic conductance produced by a presynaptic action potential as the perturbation. Separate derivations were required for modes in which the firing order is preserved and for those in which it alternates. Networks composed of two model neurons coupled by reciprocal inhibition were examined to test the predictions. The parameter regimes in which both types of nearly synchronous modes are exhibited were accurately predicted both qualitatively and quantitatively provided that the synaptic time constant is short with respect to the period and that the effect of second order resetting is considered. In contrast, PRC methods based on weak coupling could not predict 2:2 modes and did not predict the 1:1 modes with the level of accuracy achieved by the strong coupling methods. The strong coupling prediction methods provide insight into what manipulations promote near-synchrony in a two neuron network and may also have predictive value for larger networks, which can also manifest changes in firing order. We also identify a novel route by which synchrony is lost in mildly heterogenous networks.

Appendices

Appendix I: Error minimization approach to finding modes

The minimization approach for the 2:2 lockings in which the firing order does not change was as follows. The values of φ ₁₂ and φ ₂₁ were determined at each point on a grid in the (φ ₁₁, φ ₂₂) space as follows by first setting φ ₂₁ equal to P ₁{1 − φ ₁₁ + f ₁₁(φ ₁₁)}/P ₂ − f ₂₂(φ ₂₂) per the steady state version of Eq. (3) then setting φ ₁₂ equal to P ₂{1 − φ ₂₁ + f ₁₂(φ ₂₁)}/P ₁ − f ₂₁(φ ₁₁) per the steady state version of Eq. (2). In an exact solution, the steady state versions of Eqs. (1) and (4) would be satisfied, thus we utilized the quantities P ₂{1 − φ ₂₂ + f ₁₂(φ ₂₂)}/P ₁ − f ₂₁(φ ₁₂) − φ ₁₁ and P ₁{1 − φ ₁₂ + f ₁₁(φ ₁₂)}/P ₂ − f ₂₂(φ ₂₁) − φ ₂₂ as two components of error that must be minimized simultaneously. In a true solution, the value of both errors would fall to zero. All points in the (φ ₁₁, φ ₁₂, φ ₂₁, φ ₂₂) space generated using Eqs. (2) and (3) were selected if both components of the error from Eqs. (1) and (4) were below a threshold, then all adjacent points on the (φ ₁₁, φ ₂₂) grid were considered to form a cluster. The local minimum of the cluster was used as the initial condition for a gradient descent method to find the global minimum on the (φ ₁₁, φ ₂₂) grid. If the error threshold was too high, a single cluster might contain two zeroes, but the algorithm would find only one. If the error threshold is set too low, a zero might be missed entirely because none of the grid points are close enough. It is also possible for the algorithm to find a local minimum that is not zero, thus the algorithm must be applied very carefully. This algorithm also finds all 1P modes, since they also satisfy the 2P criteria in which the firing order does not change.

The minimization approach for the 2:2 lockings in which the firing order changes on every cycle was as follows. The values of φ ₁₁ and φ ₂₁ were determined at each point on a grid in the (φ ₁₂, φ ₂₂) space as follows by first setting φ ₁₁ equal to P ₂{1 − φ ₂₂    +    f ₁₂(φ ₂₂)}/P ₁ per the steady state version of Eq. (5) then setting φ ₂₁ equal to P ₁{1 − φ ₁₂   + f ₁₁(φ ₁₂)}/P ₂ per the steady state version of Eq. (7). In an exact solution, Eqs. (6) and (8) would be satisfied, thus we utilized the quantities φ ₁₂ + f ₁₁(φ ₁₁) − P ₂{1 + f ₂₂(φ ₂₁)    + f ₁₂(φ ₂₂)}/P ₁ − φ ₁₁ and φ ₂₂ + f ₁₂(φ ₂₁) − P ₁{1 + f ₂₁(φ ₁₁)    + f ₂₁(φ ₁₂)}/P ₂ − φ ₂₁ as two components of error that must be minimized simultaneously, and proceeded as described above.

Appendix II: Details of graphical method

2.1 Firing order is preserved

In order to obtain the blue curve in Fig. 4, a loop is performed over all values of ϕ ₂₂. At each value of ϕ ₂₂ chosen, the value of ϕ ₂₁ that satisfied ts₁₁ = tr₂₂, ts₁₂ = tr₂₁, and ts₂₁ = tr₁₁ was determined. The values of ϕ ₁₁ and ϕ ₁₂ were also required to determine all of the appropriate intervals. An initial estimate of ϕ ₁₁ was obtained by ignoring f ₂₁(φ ₁₂*) in the steady state version of Eq. (1), then an initial estimate of ϕ ₂₁ is made using the steady state version of Eq. (3). Then ϕ ₁₂ is estimated from the steady state version of Eq. (2). The estimate is refined by repeating the process, now considering rather than ignoring f ₂₁(φ ₁₂*) iteratively until it converges. If there was a solution, the algorithm converged in all cases tested. One problem is that it is possible that there are multiple values of ϕ ₂₁ at a given value of ϕ ₂₂, whereas this algorithm would only find one. This rarely caused a problem, however. The choice to plot tr₂₂ and ts₁₂ as shown in Fig. 4(a) is arbitrary. It is only necessary to pick two intervals for this curve that are equal to two intervals that can be calculated for the other curve as described below.

In order to obtain the red curve in Fig. 4, a loop is performed over all values of ϕ ₂₁. At each value of ϕ ₂₁ chosen, the value of ϕ ₂₂ that satisfied ts₁₁ = tr₂₂, ts₁₂ = tr₂₁, and ts₂₂ = tr₁₂ was determined. The values of ϕ ₁₁ and ϕ ₁₂ were also required to determine all of the appropriate intervals. An initial estimate of ϕ ₁₂ was obtained by ignoring f ₂₁(φ ₁₁*) in the steady state version of Eq. (2), then an initial estimate of ϕ ₂₂ is made using the steady state version of Eq. (4). Then ϕ ₁₁ is estimated from the steady state version of Eq. (1). The estimate is refined by repeating the process, now considering rather than ignoring f ₂₁(φ ₁₁*) iteratively until it converges. It should be noted that in order to interpolate properly, the intersections were actually calculated in the (ϕ ₂₁, ϕ ₂₂) plane, then the values of ϕ ₁₁ and ϕ ₁₂ were calculated at the intersection, and these values were used to predict the intervals observed in each mode. It was necessary to check that all calculated phases were in the range 0 to 1 and that the calculated intervals were nonnegative.

2.2 Firing order is not preserved

In order to obtain the red curve in Fig. 5, a loop was performed over all values of ϕ ₁₂, and the values of φ ₂₂ that satisfied ts₁₁ = tr₂₁, ts₂₁ = tr₁₁, and ts₂₂ = tr₁₂ was determined. First the value of ϕ ₂₁ was determined from the steady state version of Eq. (7). Then a loop was performed through the values of ϕ ₂₂ with ϕ ₁₁ set to the value determined by the steady state version of Eq. (5) in order to find all the values of ϕ ₂₂ that satisfied the steady state version of Eq. (8). The values of the phases so determined were used to compute tr₁₂ and ts₁₂ as shown. Multiple values of ϕ ₂₂ were sometimes found, and due to the way in which the loop is structured, points on the same branch may not be found in sequential order, but rather intermixed with points on other branches. Since the plot is a planar section of a four dimensional space, there is no guarantee that all branches are actually coplanar, thus care must be exercised in using these plots. Again, it was necessary to check that all calculated phases in the range of 0 to 1 and that the intervals are nonnegative only in a limited range of values.

In order to obtain the blue curve in Fig. 5, a loop was performed over all values of ϕ ₂₂, and the values of ϕ ₁₂ that satisfied ts₁₁ = tr₂₁ ts₁₂ = tr₂₂, and ts₂₁ = tr₁₁ was determined. First the value of ϕ ₁₁ was determined from the steady state version of Eq. (5). Then a loop was performed through the values of ϕ ₁₂ with ϕ ₂₁ set to the value determined by the steady state version of Eq. (7) in order to find all the values of ϕ ₂₂ that satisfied the steady state version of Eq. (6). The values of the phases so determined were used to compute ts₂₂ and tr₂₂ as shown. Again, the actual intersections were calculated in the (ϕ ₁₂, ϕ ₂₂) plane to facilitate interpolation.