Ramana Dodla and Charles J. Wilson  Phys Rev E Stat Nonlin Soft Matter Phys. 88:042704

Large scale simulations of electrically coupled neuronal oscillators often employ the phase coupled oscillator paradigm to understand and predict network behavior. We study the nature of the interaction between such coupled oscillators using weakly coupled oscillator theory. By employing piecewise linear approximations for phase response curves and voltage time courses, and parameterizing their shapes, we compute the interaction function for all such possible shapes and express it in terms of discrete Fourier modes. We find that reasonably good approximation is achieved with four Fourier modes that comprise of both sine and cosine terms.

Figure 1. Parameterizing the PRC and the voltage shapes.

(a) The voltage (top) and the PRC (bottom) of the HH model (thin lines) at an applied current of 10μA/cm² are taken as empirical models to construct parameterized and piecewise linear shapes (thick lines; displayed here for A′ = 0.567, B′ = −0.5, and W′ = 0.075) formulated using Eqs. 1 and 2. Parametrization of these curves allows for consideration of very general PRCs and voltage forms. (b) Comparing the interaction function and the growth function obtained from the full model (thin lines) and those obtained from piecewise linear approximations (thick lines). Fourier expansion approximations using first few modes to the H(φ) that is computed using the piecewise linear formulation are shown in dashes. The growth function, G(φ), is simply H(−φ) − H(φ), and its intersections are the phase-locked solutions of the coupled system. The negative slopes at the intersections indicate stability. The piecewise model and the Fourier approximation predicted the in-phase (φ = 0) and antiphase (φ = 1/2) states and their stability accurately. The stability of the other phase locked states is also accurately predicted.

Figure 2. Interaction function at zero spike width.

(a1) Profile of type-1 PRC with B′ = 0, at five different skewness values. (a₂) Analytically determined H(φ) (Eq. 9) for the five PRCs of (a₁). The dashed curves for A′ = 0.1, 0.3, 0.5, and 0.7 are Fourier approximations in Eq. 11. For A′ = 0.9 up to eight modes are needed, but the curve is plotted using only four modes in order to compare with the other curves. In the side panels the coefficients of cosine terms [cos(m 2π φ]), and the sine terms [sin(m 2π φ)] in the discrete Fourier spectrum of H(φ) for the first three modes are displayed for all levels of skewness. The side panel also displays FN for N = 1, 2, 3, and 4 as a function of normalized skewness. (b1) Profile of type-2 PRC at four different levels of B′, and fixed skewness of A′ = 0.3. (b2) H(φ) for the PRCs in b1 using Eq. 7 (solid) and Fourier approximations (dashed, Eq. 12). The side panels are similar to those in b2 but as a function of type parameters at A′ = 0.3. (c) The weight factor FN that captures the relative Fourier power in modes 1 to N in the parameter space of skewness and B′ as incrementally more number of modes are included in an expression for H(φ) in terms of its Fourier expansion terms. The black dots indicate the contour lines on which FN = 0.9. The brighter yellow parameter regimes bordering the dotted boundaries can be represented by the given or fewer number of Fourier modes, and contain 90% or more total power in them. For example, in the region marked F4 > 0.9 in the last panel the interaction function could be approximated by either 4 or fewer modes.

Figure 3. Interaction function at non-zero spike width.

(a1) Voltage time course with W′ = 0.075 used in all the figure panels, and profile of type-1 PRC with B′ = 0 at five different skewness values. (a2) Numerical H(φ) (solid) for PRCs in a1, and Fourier approximations for four of the curves as in Eq. 13. The approximation for A′ = 0.8 requires 5 modes but only four modes are used for comparison with the other curves. The side panels are similar to those in Fig. 2. (b1) Profiles of five type-2 PRCs at a fixed skewness of A′ = 0.3. (b2) Same as that in a2 but for the type-2 PRCs described in b1. The Fourier approximations are as in Eq. 14. (a3, b3) Fourier weight coefficient (FN) shown in the parameter space of skewness and B′ as incrementally more number of modes are included in an expression for H(φ) in terms of its Fourier modes. The black dots indicate the contour lines on which FN = 0.9 for the corresponding N. The brighter yellow parameter regimes bordering the dotted boundaries (FN > 0.9) can be represented by the given number of Fourier modes that contain 90% or more total power in them. Type-1 PRC is illustrated with B′ = 0 in a3 and type-2 PRC with B′ = −0.5 in b3.

Figure 4. Proportion of odd components in H(φ).

Oddness factor F_odd in the plane of skewness (A′ that ranges from 0 to its maximum allowed value of 1 − W′) and type parameter (B′ shown here between −1 and 1) at four different values of spike width. The interaction function is rarely totally even or totally odd. Even predominantly evenness (F_odd < 0.1) or predominantly oddness (F_odd > 0.9) occurs in tiny pockets of the parameter space. For most of the type parameter and skewness values H(φ) is a mixture of both even and odd modes.