A Model of the Effects of Applied Electric Fields on Neuronal Synchronization

2.1. Two-Compartment Neurons Embedded in a Resistive Array

We based our numerical experiments on a computational model that we have explicitly designed to include electric field interactions (Gluckman, 1998). The presence of an electric field induces spatial polarization in neurons (Chan and Nicholson, 1986; Chan et al., 1988; Tranchina and Nicholson, 1986), and therefore, the minimal individual neuronal unit must have at least two spatially separated compartments. We therefore chose the two-compartment Pinsky-Rinzel (PR) model neuron, which consists of a dendrite and a soma compartment separated by a finite conductance g_c (Pinsky and Rinzel, 1994). A schematic representation of the PR model neuron is shown in Fig. 1(a). The dendritic compartment contains a passive leakage current, a capacitive current, and three voltage-dependent ionic currents. The soma compartment contains only two voltage-dependent ionic currents, as well as passive leakage and capacitive currents. Steady currents can also be injected independently into each compartment.

For the present work, two such neurons are embedded in parallel within an array of resistors which models the extracellular medium; see Fig. 1(b). The terminal resistors on both ends of the array are chosen to match the impedance characteristics of an array of infinite extent. An externally applied electric field is modeled by imposing a potential difference V_app between point ‘A’ and ground. The existence of this resistive array also allows the neurons to communicate ephaptically. Finally, the neurons are synaptically coupled to each other via both NMDA and AMPA channels. Complete details regarding the resistive array, various currents, rate functions, and other model parameters are presented in the Appendices.

The focus of our investigation is to study the synchronization properties of a heterogeneous network of neurons under the influence of an applied electric field. To accomplish this goal, we modulate the applied potential difference V_app and we adjust the degree of heterogeneity within the network by varying the internal conductance g_c between the neurons’ somatic and dendritic compartments. The entire system of differential equations was integrated using a 4th order Runge-Kutta algorithm with a fixed integration time step.

2.2. The Measurement of Phase from Spike Times and the Synchrony Measure

Throughout this study, we specifically define neuronal synchrony in terms of the phase locking between the neurons’ spiking activity. In general, uncoupled heterogeneous neurons will spike at different rates. When these neurons are coupled together, one would expect these different spiking frequencies to “pull” toward a common mean frequency if the heterogeneity is not too large (Winfree, 1980; Kuramoto, 1984). The spike timings of the coupled neurons are expected to “phase-lock” to each other with a possible non-zero phase lag Ψ₀. With still stronger coupling, the (non-identical) neurons might also spike at nearly the same time, i.e., Ψ₀ → 0, in addition to spiking at the same rate. In this study, we consider two heterogeneous neurons to be synchronized if they simply phase lock to each other. We do not impose the stronger condition of exact synchrony with zero phase lag. The goal of this study is to examine the relationship between the phase-locked state, the applied electric field, and the degree of heterogeneity among the neurons.

Our procedure for assigning a phase to a given neuron’s activity is as follows. A voltage threshold for the somatic voltage V_s is chosen, and the spike times t_n , where n = 1, 2, 3. . . , are defined as the times when V_s crosses this threshold.¹ Then, the phase corresponding to the activity of neuron i at an arbitrary time t between its nth and (n + 1)th spike (t_n⁽ⁱ⁾ ≤ t < t_n+1⁽ⁱ⁾) is defined as ${\phi }_i(t)=2\pi (\frac{t-t_n^{(i)}}{t_{n+1}^{(i)}-t_n^{(i)}})+2\pi n$. Thus, the phase increases linearly by 2π between successive spikes (see Fig. 2(a) and (b)).

A general n:m phase-locked state can be defined by the condition |nφ₂ − mφ₁ − Ψ₀| < ɛ, where n and m are integers (n, m = 1, 2, 3 . . .), φ₁ and φ₂ are the phases of the two neurons, Ψ₀ is a constant phase lag between 0 and 2π , and ɛ is a small constant (Rosenblum et al., 1996; Pikovsky et al., 2000).² For the parameter range that we examined, the 1:1 phase-locked state was observed predominately, i.e., |φ₂(t) − φ₁(t) − Ψ₀| < ɛ. In our analysis, we introduce the relative phase Ψ(t) = φ₂(t ) − φ₁(t) between the two neurons. From the definition of phase given above, Ψ(t) is simply the instantaneous phase of neuron-2, φ₂(t ), when neuron-1 completes its k-th cycle at time t = t_k⁽¹⁾:

In terms of this relative phase, the degree of phase-locking can be quantified by the synchronization indexn γ (Kuramoto, 1984; Pikovsky et al., 2001), defined as

where 〈〉 is an average over the number of spiking events (N ), which we set to 10000 in most of our simulations. If one assigns a two-dimensional unit vector (a phasor) to each relative phase measurement of Ψ, then γ gives the magnitude of the vector sum of all the phasors (see Fig. 2(c)). One can imagine that in the incoherent case, all measurements of Ψ(t_k⁽¹⁾) and the corresponding phasors will be uniformly distributed on a unit circle and γ should approach zero for large N. On the other hand, if the neurons are phase-locked, the distribution of the relative phase Ψ(t_k⁽¹⁾) will cluster around the locked value, Ψ₀, and γ should approach its maximal value of one (see Fig. 2(c)).

3.1. Natural Frequency Mismatch

In this section, we present results on the synchronization characteristics of our coupled neurons as the strength of the externally applied electric field and the heterogeneity parameter Δg_c are varied. Figure 3 illustrates the variation of the synchrony index γ as the applied electric field is modulated in the case of two nearly identical neurons (solid squares) and two very different neurons (open circles). For inhibitory (negative) electric fields, the neurons tend to phase-lock with each other, and there is a range of moderate inhibitory field where the spiking ceases (see gap in the figure). Larger negative fields can depolarize the dendrites and paradoxically “excite” the neurons back into phase-locked spiking behavior.³ Qualitatively, this response to negative inhibitory fields is consistent across different levels of parameter mismatch.

However, with excitatory (positive) electric fields, we observed a more complex network synchrony response as a function of the applied electric field and of heterogeneity. Neurons with a small parameter mismatch first desynchronize at moderate excitatory fields, and then resynchronize at larger excitatory fields.⁴ Neurons with a large parameter mismatch become desynchronized as the strength of excitation becomes stronger, and remains so. This result is consistent with the report of Golomb and Hansel (2000), who found that in sufficiently heterogeneous ensembles, increasing coupling tends to desynchronize the network. We show in Figs. 3(b) and (c) tracings of the somatic voltage from the two neurons at large excitatory fields (E = 600 mV/cm) in the synchronous and the asynchronous cases just described. In the synchronous case (Fig. 3(b)), the spiking from the two neurons phase-lock, and this behavior is reflected in the constant relative-phase (solid dots) plotted in the graph. In Fig. 3(c), we illustrate the asynchronous case in which the relative phase between the two neurons monotonically drifts in time.

A more complete illustration of our network’s synchronous behavior is presented in Fig. 4(a). In this figure, we show the γ index as a function of both the applied electric field strength and the parameter mismatch Δg_c. We discuss several prominent features. An extensive asynchronous region is present for large excitatory fields and large neuronal disparity, in agreement with Golomb and Hansel (2000). Moderate inhibitory fields tend to induce phase locking even when neuronal disparity is large. Larger inhibitory fields suppress collective network activity. This is consistent with in vitro seizure suppression experiments (Gluckman et al., 1996, 2001). Lastly, the reemergence of phase-locked spiking activity under the largest negative fields is also seen. Among these general characteristics, the most interesting feature of this phase diagram is the boundary between the phase-locked and asynchronous (phase-drifting) states under an applied excitatory field. For small neuronal disparity there is a regime (see the wedged region indicated in Fig. 4) where moderately strong excitatory fields desynchronize the network, but an even stronger field leads to resynchronization. The threshold for this transition at small neuronal disparities is seen in Fig. 4(a) to increase with increasing Δg_c.

We now argue that the transition to synchrony in the network including the wedged region can be explained using a simple phase oscillator model (Winfree, 1980; Cohen et al., 1982; Kuramoto, 1984; Ermentrout and Kopell, 1984). Specifically, the boundary between the synchronous and the asynchronous states is highly correlated with the degree of natural frequency mismatch among the individual units within the network. To illustrate this connection, we plotted the degree of natural frequency mismatch Δω between the neurons in isolation as a function of the applied electric field and neuronal disparity in Fig. 4(b). For each neuron in isolation, its intrinsic natural frequencyω is measured with the same parameters as in the coupled network and subjected to the same applied electric field. The natural frequency mismatch Δ ω is then simply the difference in theω values between the neurons. By comparing Figs. 4(a) and (b), one can see that there exists a critical value of natural frequency mismatch Δω^* such that for Δω > Δω^*, the network is in the asynchronous state, and for Δω < Δω^*, the neurons phase lock to each other. The critical value Δω^* was empirically determined in our numerical experiment to be approximately 0.63 rad/s. It is remarkable that this simple criterion reproduces the observed synchrony transition boundary so well. In particular, it correctly predicts that for small neuronal disparity (Δg_c < 6%), increasing the (excitatory) electric field first desynchronizes and then resynchronizes the network (see the wedged region in Fig. 4(b)). To clarify this point, Fig. 5 (lower panel) shows the decrease in Δω as a function of the increasing (excitatory) electric field strength at Δg_c = 4%. One can understand this result by noting that each individual neuron’s spiking rate is a monotonically increasing concave-down function with respect to the applied electric field (see Fig. 5—top panel). This behavior is biologically reasonable since the spiking rate of a neuron will ultimately be limited. Thus, as the spiking rate increases with increasing excitatory fields, the corresponding incremental change between two different neurons is smaller. In the small disparity regime, an initial frequency mismatch might be larger than Δω^* at low excitatory fields (so that the network is asynchronous when coupled), but it might become small enough with increasing excitatory fields such that the difference between the spiking rates of the two neurons decreases below the critical level Δω^*. Thus, at sufficiently high excitatory field, phase locking is observed.

In our discussion so far, we have focused on the simplest 1:1 phase locked state. Within the parameter range of our study, we have observed a few cases where the two neurons are synchronized at higher-order n : m phase locked states (cross-hatched squares in Fig. 4(a)). If Ω₁ and Ω₂ are the observed frequencies of the coupled neurons, the n : m phase locked state is characterized by the condition: nΩ₂ − mΩ₁ = 0 (Pikovsky et al., 2001). These higher order locked states are expected to occur when the natural frequencies of the two neurons in isolation are roughly in the same rational ratios, i.e.,ω₁/ω₂ ≈ n/m (Ermentrout, 1981; Pikovsky et al., 2001). Figure 6 gives some examples of the values of Δg_c (10%, 14%, and 24%) at E = −100 (mV/cm) where the natural frequencies of the neurons give the approximate rational ratios 3/4, 2/3 and 1/2. To demonstrate these n:m phase locked states, we plotted ΔΩ = nΩ₂ − mΩ₁ as a function of Δg_c for four different combinations of n and m (see Fig. 7). Since Δg_c is functionally related to Δω and is a direct measure of the difference between the individual uncoupled neurons, Δg_c serves as the detuning parameter for this analysis. The n:m phase locked state in these detuning curves are indicated as plateau regions with Δ Ω = 0. As one can see from Fig. 7, the plateau regions correspond to ranges in Δg_c where ω₁/ω₂ ≈ 1(0% ≤ Δg_c ≤ 4%), 1/2 (22% ≤ Δg_c ≤ 26%), 2/3(13% ≤ Δg_c ≤ 14.2%), and 3/4 (9.9% ≤ Δg_c ≤ 10.3%).⁵

3.2. Phase Response

In our network, the dynamics of the individual units are determined by the Pinsky-Rinzel equations, which include many state variables. Nevertheless, the resulting dynamics in our chosen parameter range are predominately periodic spikes. This periodic behavior can typically be reduced to a simple phase equation under a suitable parameterization and the assumption of weak coupling. Various authors (Rand and Holmes, 1980; Cohen et al., 1982; Kuramoto, 1984; Ermentrout and Kopell, 1984) have given detailed descriptions for similar phase reductions in different biological systems. For the following discussion, we assume that the reduced phase equation for the ith uncoupled neuron can be written as

where φ_i is the phase of the ith neuron, and ω_i (E, g_c⁽ⁱ⁾) is the natural frequency of the ith neuron in isolation (i.e., without any synaptic or electrical couplings between neurons). Note that ω_i depends explicitly on the externally applied electric field E and the internal conductance g_c.

With two or more neurons mutually coupled together as in our network, we follow Refs. (Kuramoto, 1984; Cohen et al., 1982; Ermentrout and Kopell, 1984) and further assume that the input from the presynaptic neuron is not so strong that it significantly alters the form of the spike. From our numerical simulations, this is a valid assumption for almost the entire range of parameters examined. Under this assumption, the network of coupled phase oscillators can be written as

where the sum is over all neurons (j) which are connected to the ith neuron. In general, we expect the coupling function f_i (E, φ_j − φ_i) to depend on the applied electric field E and the relative phase of the neurons. In the case of only two mutually coupled neurons, we can explicitly write down the following reduced dynamical equations:

where Ψ = φ₂ − φ₁ is the relative phase between neuron-1 and neuron-2.

To understand the dynamics of phase locking within this network, one can consider the temporal evolution of the relative phase Ψ by subtracting the previous two equations,

where Δω ≡ ω₂ − ω₁ is the mismatch in the natural frequencies of the two individual neurons in isolation and Δf ≡ f₁ − f₂, the phase sensitivity function, is simply the difference between the coupling functions. For identical neurons with symmetric coupling, f₁(ψ) = f₂(− ψ) = f (ψ) and Δf (ψ) will simply be twice the odd part of f (ψ), i.e., 2 f_odd(ψ) = f (ψ) −f (−ψ). For non-identical neurons, Δf (ψ) will typically contain both even and odd parts. For a pair of neurons to phase lock to each other, the necessary condition is $\dot{\mathrm{\Psi }}=0$, i.e., the relative phase Ψ is constant. For a particular choice of the parameters E and Δg_c, the natural frequency mismatch of the neurons, Δω, is constant and Δf is a function of the phase difference Ψ only. Requiring the right hand side of Eq. (2) to be zero, the condition for phase-locking becomes

where Ψ* is the locked phase lag between the two neurons. A schematic illustration of this phase-locked criterion is given in Fig. 8. The dashed line is the constant value of Δω and the solid line is the function Δ f (Ψ). The crossing points (indicated by the three small circles) are the possible phase-locked equilibrium states predicted by this reduced phase model for a particular choice of applied electric field E and Δg_c. The stability of these phase-locked states is given by the sign of the local slope of Δf (Ψ) at these locations. For example, Ψ*_b is a stable equilibrium with $\frac{d\mathrm{\Delta }f(\mathrm{\Psi })}{d\mathrm{\Psi }}{|}_{{\mathrm{\Psi }}^{*}}>0$. One can understand the stability argument by noticing that for Ψ slightly less than ψ*_b^, Δf (ψ) < Δω so that $\dot{\mathrm{\Psi }}>0$. Then, in time, Ψ will increase toward Ψ*_b from the left. Similarly, for Ψ slightly larger than Ψ*_b, Δf (Ψ) > Δω and $\dot{\mathrm{\Psi }}<0$ so that Ψ will decrease toward Ψ*_b from the right. By a similar argument, one can see that the remaining two phase-locked states at ψ*_a and ψ*_c are unstable equilibria. In summary, if one can obtain the phase sensitivity function Δf (Ψ) and the natural frequency mismatch Δω for the ensemble at a given choice of E and Δg_c, one can then use the reduced phase model to predict the stable phase lags Ψ*_b for the two locked neurons from a graph similar to Fig. 8.

To get a conceptual picture of the synchrony transition for a pair of heterogeneous neurons, we discuss a canonical example for illustrative purposes. Consider the simple example in which f and Δf are given by Figs. 9(a) and (b). The coupling function f is basically given by its canonical form of (a₀ + a₁ cos(ψ) + b₁ sin(ψ)) (Ermentrout, 1996) and Δf is simply the anti-symmetric part of f. Furthermore, since the elicitation of spikes in one neuron from another neuron cannot be instantaneous, this coupling function f will typically be phase shifted by a small amount δ. One should in general expect this phase shift to be influenced by the synaptic time constants of the NMDA or the AMPA channels, the time constant in the kinetics of the membrane, and the conduction velocity between neurons (not modeled here).

To analyze the phase-locked states of this idealized model, we need to compare Δf (Ψ) with Δω. For identical neurons, the natural frequency mismatch Δω is zero. Then according to Eq. (3), there will be two phase-locked equilibria at Ψ* = 0(2π ) and π (Fig. 9(b)). Of these two phase-locked states, only the solution at Ψ^* = π has a positive slope. Thus, Ψ^* = π (the anti-locked state) is the only stable equilibrium. For non-identical neurons, Δω will be non-zero and its value will depend on the parameter mismatch Δg_c between the two neurons and the applied electric field E. Depending on the sign of Δω, the relative phase Ψ^* of the stable phase-locked state will either decrease (Δω < 0) or increase (Δω > 0) away from the anti-locked phase. There are two special relative phase values Ψ_L^c and Ψ_R^c corresponding to the left (negative) and the right (positive) peaks of the Δf (Ψ) function. These are, respectively, the smallest and the largest locked relative phases possible for the system. The critical values for this coupling threshold are achieved when Δω = Δf (Ψ_L^c) or Δf (Ψ_R^c). These two values also define the range of allowed frequency mismatch (Δf (Ψ_L^c) ≤ Δω ≤ Δf (Ψ_R^c)) where phase-locking between the two neurons is possible. For Δω outside this range, a non-locked phase drifting state between the two neurons is expected. One should note that, for non-identical neurons, an exact phase-locked state with Ψ = 0 is not expected to be stable in this simplified model. Biologically, this is a reasonable conclusion since synaptic coupling cannot be instantaneous, and no two neurons are exactly identical.

The key advantage of this phase model description is that it predicts the synchronization characteristics of the coupled network from the knowledge of the individual neurons, namely: their natural frequencies and coupling functions. We now apply this formalism to our model system. For a given level of applied electric field E and internal conductance g_c, the natural frequency ω(E, g_c) for a particular neuron can be determined simply by taking the inverse of the average of the spiking rate.⁶ The coupling functions f_1,2(E , φ_j − φ_i) in the reduced phase model can be rigorously calculated using the averaging technique described in Refs. (Ermentrout and Kopell, 1984; Hansel et al., 1995; Ermentrout, 1996). In brief, suppose a pair of neurons is described by the following dynamical equations:

where X is a multi-dimensional vector describing the state of the neuron, F₁ and F₂ prescribe the intrinsic dynamics of the neurons in isolation, and G₁ and G₂ are the coupling functions between the two neurons. It can be shown that when the long-term dynamics of the individual neurons are stable limit cycles and the coupling strength ɛ is small, then the above nonlinear equations can be reduced to a pair of phase equations as in Eq. (1). The coupling functions f_1,2 in the phase equations can then be expressed as a time-averaged quantity involving the original coupling functions G_1,2 as follows:

where X₀(t ) is the period-P limit cycle of neuron j in isolation and X^*(t), called the adjoint solution, is determined by the linearized dynamical equation

where [D_X F]^T is the transpose of the Jacobian matrix of F(X) with respect to the state variable X. X^*(t) is then further normalized according to the condition X*(t) · X′₀ (t) = 1, where the prime denotes the rate of change of the vector field along the periodic orbit X₀(t). This procedure is elegantly implemented in the publicly available software package XPPAUT (Ermentrout, 2002).

Alternatively, one can measure the coupling function f (E, φ_j − φ_i) of one neuron with respect to another by a sampling method. One can employ a pair of unidirectionally coupled subsystems. For two mutually coupled neurons, we decompose the system as shown in Fig. 10. In Fig. 10(a), neuron-1 receives synaptic input from neuron-2 only and one measures the rate of change of neuron-1’s phase, ${\dot{\phi }}_1$ concurrently with φ₁ and φ₂. Then, the quantity ${\dot{\phi }}_1-{\omega }_1$ gives an estimate of f₁(E, φ₂ − φ₁) for a particular phase difference Ψ = φ₂ − φ₁. In some cases, the two neurons phase lock with each other under uni-directional coupling. In this latter situation, one can still sample f₁(E, φ₂ − φ₁) by periodically perturbing the pair of neurons by a sufficiently large applied electric field or by injecting current, and measuring ${\dot{\phi }}_1$ during the system’s transient back toward the phase-locked state.

Figure 11 is an example of the coupling function calculated using Eq. (4) for a PR neuron embedded within a resistive network with an excitatory applied electric field. The inset is a close approximation to this function using the sampling method described above. Both of these methods agree qualitatively with the canonical functional form of the coupling function, and the δ phase shift due to synaptic delay can be clearly seen in the graph.

Since the neurons are coupled both synaptically and electrically through the resistive network, one can in principle tease out their separate contributions to f (ψ). This can easily be done by utilizing Eq. (4) and taking the time average of, respectively, the synaptic current I_syn and/or the effective current term g_c V_DS^out resulting from the electric polarization between the two chambers (see Appendixes A and B for details on these two terms). In Fig. 12, we have plotted the coupling function for a neuron with the same parameters as in Fig. 11 but with no externally applied electric field. In this case, we have chosen no applied electric field because we simply want to examine the phase response of a neuron purely due to the firing of its neighbor within the resistive array. Figures 12(a) and (b) are, respectively, the graphs for the coupling function f (ψ) with and without I_syn. Since the neuron remains embedded within the resistive array, the firing of the neighboring neuron can still have a modulating effect on its firing time even when I_syn = 0. However, as one can see from this graph, the contribution to the coupling function from the ephaptic coupling (Fig. 12(b)) is substantially smaller than the contribution from the synaptic channels (Fig. 12(a)). Nonetheless, since it is the difference between the coupling functions Δf (E, Ψ) that is important in determining the phase-locked criterion, this small contribution plays a significant role in determining the synchronization properties of the neurons. Furthermore, it is important to note that in contrast to the larger but positive-definite contribution from the excitatory synapes, the ephaptic contribution can be both positive and negative. Thus, the ephaptic connection has the potential to both advance and retard the phase of the next spike of the postsynaptic neuron.

To illustrate the utility of the reduced phase model in analyzing the synchronization among our PR neurons, we consider three examples respectively in the locked in-phase state, the phase drifting state, and the locked anti-phase state.

Our example for the in-phase locked state is from two non-identical phase-locked neurons. The parameters were E = 600 mV/cm, g_c⁽¹⁾ = 2.1 mS/cm², and g_c⁽²⁾ = 2.058 mS/cm². From the direct simulation of this system, the ensemble has a single locked state with a phase-lag between the two neurons near Ψ = 5.7 rad (see Fig. 13(c)). The coupling functions f_{1, 2} for both neurons were calculated using XPPAUT and the results are graphed in Fig. 13(a). The phase sensitivity curve Δf for these two coupled neuron is shown in Fig. 13(b). Note that it crosses the natural frequency mismatch Δω = 0.5 rad/sec line near its right positive peak. The stable equilibrium state predicted from this crossing point roughly agrees with the phase locked value of Ψ = 5.7 rad observed in the full model (Fig. 13(c)). As we have discussed previously, as Δω increases, the stable equilibrium moves toward the in-phase regime near Ψ = 0 (or 2π ). Due to the intrinsic synaptic delays, the locked phase between the two neurons cannot be exactly at Ψ = 0 (or 2π). This slight phase lag in the in-phase locked neurons is exactly what we have observed in our simulation.

The next example is for the case when the two neurons do not phase lock to each other. The parameters were chosen in the asynchronous region with E = 200 mV/cm, g_c⁽¹⁾ = 2.1 mS/cm², and g_c⁽²⁾ = 1.68 mS/cm². In this case, the phases of the two coupled PR neurons drift with respect to each other, as can be seen in Fig. 14(b). Again, using XPPAUT, we have calculated the coupling functions for these two neurons in the reduced phase model and the phase sensitivity graph Δf is plotted in Fig. 14(a). As one can see from this graph, Δω does not intersect Δf at all, and the two neurons are not expected to phase lock. This is what is observed in the actual numerical experiment.

The last example is for the anti-phase locked state. As we have seen in our simplified phase model (Fig. 8) previously, identical neurons (Δg_c = 0%) with excitatory coupling typically phase lock in the anti-phase state (Hansel, 1995). For E = 200 mV/cm and g_c^(1,2) = 2.1 mS/cm², Δω is zero and Δf is given in Fig. 15(a). As expected, this curve crosses Δω = 0 near Ψ = π. However, upon magnification, one can see that Δf actually has three zeros in this region with the left and right zeros being stable and the middle one at ψ = π being unstable. Therefore, the reduced phase model predicts that there should be two distinct locking regions near the anti-locked state. We ran many trials (2000) of our numerical experiment using the full model of the coupled PR neurons and formed a histogram of the observed locked relative phases (see Fig. 15(c)). From this histogram, we found two main clusters on either sides of Ψ = π which indicate a multi-stable phase-locked state in this anti-phase regime. The two peaks of the histogram align closely with the two predicted locked phases. However, if one examines the histogram closely, one notices that there is a spread to the two histogram peaks. This indicates additional structure in the anti-synchronous regime for the two-coupled PR neurons. We speculate that the ephaptic coupling through the intercellular medium may play a subtle role in modulating the synchronization between these two neurons within this multi-stable region.