neuronal synchrony is implicated in working memory, sensory processing, and attention (Fries et al. 2001; Hansel and Sompolinsky 1996; Wang 2001) and diseases, such as epilepsy, Parkinson's, schizophrenia, Alzheimer's, and autism (Uhlhaas and Singer 2006). Cortex and hippocampus, due to their laminar structure, which results in large local field potentials (LFPs), are major research objects in the study of neuronal synchrony and epilepsy (Dudek and Sutula 2007; Jiruska et al. 2010b; McCormick and Contreras 2001). Neuronal network hyperactivity and hypersynchrony underlie activity bursts and spreading excitation, which are related to transient seizure episodes (interictal events) and seizure onset, respectively (Jefferys 1995; Jiruska et al. 2010b; Traub et al. 2004; Weiss and Faber 2010).

The most widely discussed mechanism for pathologic hypersynchrony in neuronal networks is the remodeling of synaptic connections. Experiments show that loss of inhibitory interneurons and the formation of recurrent connections can lead to hypersynchrony (Dudek and Sutula 2007). Whereas synaptic remodeling is likely involved in epileptiform behavior, hypersynchrony can also exist without synaptic connections (Jefferys and Haas 1982). For example, extracellular medium osmolality changes induce or suppress epileptiform behavior in the absence of functioning synapses (Bikson et al. 2003; Dudek et al. 1990; Haas and Jefferys 1984; Jefferys 1995). Even though theoretical and ultrastructural studies have suggested that gap junctions between neurons could explain these observations (Hamzei-Sichani et al. 2007; Traub et al. 2000, 2004), synchrony can even persist when gap junctions are functionally blocked (Jiruska et al. 2010a). Therefore, nonsynaptic, not gap-junction-mediated, mechanisms for neuronal hypersynchrony likely exist.

One potential mechanism is ephaptic coupling, meaning coupling via endogenous electric extracellular potentials (V e ). Low-frequency V e , simulating spatially averaged potentials from many cells, i.e., LFPs, can synchronize neurons (Anastassiou et al. 2011; Berzhanskaya et al. 2013; Frohlich and McCormick 2010; Radman et al. 2007). It is unclear, however, how similar-magnitude V e , generated by single action potentials, affect population spike timing. Because V e generated by individual spikes decay quickly with distance and are highly transient, single-spike ephaptic coupling is often disregarded (Buzsáki et al. 2012). However, a theoretical study predicted that single-spike ephaptic coupling might significantly alter spike times between a pair of neurons (Park et al. 2005). A single action potential can generate both hyperpolarizing and depolarizing V e , dependent on location (Gold et al. 2006, 2007). Ephaptically induced changes to the membrane potentials of nearby neurons show a similar biphasic shape (Vigmond et al. 1997). It remains unclear how ephaptic potentials and their biphasic signature affect spike timing and network synchronization (Jefferys 1995; Weiss and Faber 2010).

Here, we investigate the effect of single-spike ephaptic coupling by first modeling ephaptic coupling between neuron pairs, demonstrating its effects in a small network, and then, generalizing pairwise coupling characteristics to an oscillator network. Our results provide a proof of principle that single-spike ephaptic coupling synchronizes neurons dependent on neuronal spacing and ephaptic potential shape and underlies formation of locally synchronous microclusters.

where Δis the smallest time difference from the time of spikein spike trainto any spike in spike train, andis the number of spikes detected in spike train. The number of clusters was chosen by maximal modularity using the similarity measure

To adapt the oscillator network model to experiments with hippocampal slices, an elongated geometry of the oscillator field was introduced ( Jiruska et al. 2010a ; Muldoon et al. 2013 ). The square area was changed to a rectangular area whose long side is 10 times the length of its short side, while preserving the relationship that oscillator density was equal to 1/ R 2 (e.g., see Fig. 9 A ).

Simulations were executed by iteration through events (transition points in piecewise linear potentials, ϕ n = 1 threshold crossings that reset ϕ n to 0 and end of refractory period), increasing numerical efficiency and inherently eliminating step-size problems of finite time-step methods (e.g., Euler or Runge-Kutta schemes). It is important to realize that different from the prior analysis of conditions for synchrony, this oscillator network model was based on continuous time and the detailed interaction of ephaptically coupled oscillators throughout the entire phase cycle. This treatment was no longer based on changes in phase difference after a whole-phase cycle of the stimulating neuron.

is the time at which the ephaptic potential deflection first returns to 0. The coefficientexpresses the amplitude of ephaptic coupling between oscillators at a distance of 1 μm, which was derived from our NEURON results. It was assumed that= 0 corresponds to the resting potential and= 1 to the firing threshold for an action potential. For thepyramidal neuron simulated in NEURON, this corresponds to an ≈20-mV difference in membrane potential ( Table 1 ). Therefore, the amplitude of the change ininduced by an ephaptic potential at 1 μm was/20 mV.describes the peak-to-peak deflection ofwhen receiving and stimulating neurons are offset by 1 μm, as calculated in NEURON. To derive the ephaptically induced rate of phase change, this amplitude was divided by the time necessary to reach the ephaptically induced deviation,/2. Thus

with= 1, . . . ,. The intrinsic frequencies () were randomly drawn from a normal distribution centered onwith an SD of. The other terms represent ephaptic coupling of other neuronsto neuron, which were summed to give the overall alteration ofby ephaptic coupling.is the Euclidean distance between neuronsand. The polarity of the ephaptic-like coupling was set by= ±1 (i.e.,= +1 for “up-down” or= −1 for “down-up”;= 0 was used whento prevent oscillators from coupling to themselves).is analogous to parametersand̄ in that all determine the polarity of ephaptic coupling, althoughdoes not control the amplitude of coupling.was assigned randomly for each pairwise connection, with a probabilitythat= +1 (e.g.,= 1 means that all ephaptic potentials have an up-down shape).) is a Heaviside function, ensuring that the oscillatorwas not affected by ephaptic potentials during a refractory period of durationafter crossing= 1. Outgoing ephaptic potentials were triggered when= 1 was crossed. At the same time,was reset to zero. The three parameters, andwere changed among our different simulations. The fixed simulation parameters are shown in Table 1

Let p ( δ ) be the distribution of the phase differences of N → ∞ oscillators relative to a reference phase ϕ 0 (all oscillators and the reference phase have the same intrinsic phase frequency). To add a randomizing influence to the ephaptic-like coupling effects, a Fokker-Planck formalism is used ( Risken 1984 )

for allis the condition for a global steady state. For any two oscillators (and) from the network,must hold to fulfill the condition from. With the introduction of a third oscillator (), the possibility of a network-wide steady state can be disproved by the following contradiction with

For a network of N identical oscillators, now with bidirectional, ephaptic-like coupling between all oscillator pairs, the condition for a network steady state is that for all a ≠ b ( a , b ∈ 1, . . . , N }), the δ a,b take on a steady-state value of δ a,b * simultaneously. Based on the fraction of up-down ( c > 0) potentials ( p + ), two extreme cases can be investigated analytically. For p + = 1, all c > 0, and therefore, all δ a,b = δ a , b * = 0; all phase differences vanish, and the network-stable state is perfectly synchronous. For the opposite case, p + = 0, all c < 0, so

With the consideration of the relationship of a and b , the horizontal axis (see Fig. 2 E ) corresponds to − δ in our formulation, and Δ t AP on the vertical axis can be seen as roughly proportional to − f . At a specific zero crossing point ( ϕ e ), the biphasic f changes from positive to negative (for up-down potentials). For δ < ϕ e , the ephaptic potential from a leads to a predated firing of b (Δ t AP < 0), which corresponds to a reduction of δ ( f > 0 for up-down potentials). For δ > ϕ e , the ephaptic coupling leads to an increase of δ ( f < 0 for up-down potentials). With the assumption of a continuous and smooth f that is zero only at δ = 0 and δ = ϕ e , this means that for up-down potentials, δ* = 0 is a unique, stable steady state, attracting all initial values except δ * = ϕ e , the location of a unique, unstable steady state. With the maintenance of the same f , setting c < 0 captures down-up potentials, δ* = 0 is unstable, and δ* = ϕ e is stable. In summary, the globally attractive, stable steady states are

We constructed an analytical mathematical model consisting of a stimulating neuron ( a ) and a receiving neuron ( b ). Both neurons are modeled as simple phase oscillators with phases ϕ a and ϕ b , respectively. Neuron a fires tonically with a constant frequency, so that ϕ a is constant. ϕ b increases with the same intrinsic frequency but can be altered due to ephaptic potentials received from a , ϕ ̇ b = ϕ ̇ a + cf . Therefore, δ a,b undergoes the following time evolution

We calculated Pearson correlation coefficients between somatic V m to quantify synchronization between pairs of neurons in the NEURON network ( Poulet and Petersen 2008 ). Neurons were simulated for 500 ms, and V m from the last half of simulations were used to compute correlation coefficients. Neurons were initialized with different V m , which was assigned randomly between −70 and −40 mV. All sections in a given neuron were initialized to the same V m . All other initial conditions were identical between neurons.

We additionally implemented ephaptically coupled networks inside NEURON. Ephaptic stimulation was achieved identically to the above pairwise stimulation. V e , at each section of each neuron, was calculated by summing the contributions to electric potential from all neurons except the neuron being stimulated. Networks consisted of 16 identical layer V pyramidal cells arranged in a square grid. Spacing was varied to control the strength of ephaptic coupling. Since we used neuron spacings that were less than the extent of dendritic arbors, neuron segments could possibly overlap, although in practice, a minority of segments did (<1/2,000 of segment-to-segment pairings). Overlapping segments were not included when calculating V e .

The V e generated by one (“stimulating”) neuron was applied to a second (“receiving”) neuron. This was achieved by setting the extracellular potential of the receiving neuron equal to V e using the extracellular mechanism in NEURON ( Carnevale and Hines 2006 ). A unique V e was calculated for each section of the receiving neuron, and by calculating V e at the appropriate locations x , we could vary the proximity and relative position of the two neurons. Except for modifying V e , the receiving neuron was implemented in NEURON identically to the stimulating neuron. Without ephaptic stimulation, the receiving neuron generated a spike identically to the stimulating neuron. In the pairwise interactions, the electric potential due to the receiving neuron was not accounted for. Its effect is predicted to be small, as the receiving neuron fires during the refractory period of the stimulating neuron. We neglected the effect of a neuron's V e on itself ( Holt and Koch 1999 ).

whereis the potential due to somatic currents (). Since the soma is assumed to be spherical,can be calculated aswhereis the distance betweenand the center of the soma. Each model section was approximated as composed of multiple straight cylinders, and total current from the section was uniformly distributed across all cylinders in the section. Individual potentials () were calculated for each cylinder. For full details, see Holt (2008)

We implemented a reconstructed layer V pyramidal cell model with 176 sections in the NEURON environment using procedures and parameters of Mainen and Sejnowski (1996) . This is a multicompartment model with a high density of Na + channels in the axon initial segment and hillock. The soma and dendrites contained slow K + channels, and fast K + channels were present in the axon and soma. Spontaneous spiking was achieved by increasing synaptic conductance by 0.2× the leak conductance and setting the reversal potential to 0 mV in the dendrite compartments ( Holt and Koch 1999 ). To achieve repetitive firing without slower dynamics, Ca 2+ currents were removed from all compartments. The resulting membrane currents caused the cell to fire with a regular interspike interval of ∼15 ms. All membrane currents (capacitive and ionic) were calculated every 0.01 ms. The bulk extracellular resistivity ( ρ e ) was 400 Ωcm.

RESULTS

Biphasic V e with Millivolt Amplitude Around a Spiking Neuron We used a previously published NEURON model of a layer V pyramidal neuron to investigate ephaptic coupling from a stimulating to a receiving neuron (Mainen and Sejnowski 1996) (Fig. 1A). As per Holt and Koch (1999), the model was modified to produce an action potential. We calculated the V e , resulting from membrane currents in the stimulating neuron. The action potential began at the soma hillock and axon initial segment, sections of the neuron with Na+ channel densities 100-1,000 times greater than other areas, and propagated to the dendrites. Fig. 1.Characterizing electric potentials around a spiking neuron. A: multicompartment, 3-dimensional, layer V pyramidal cell model (Mainen and Sejnowski 1996). Duplicate model cell is shown in gray. B: summed membrane current for all compartments of the apical trunk, soma [including the hillock and initial axon segment (i.seg)], and the basal dendrites. C: 2 characteristic extracellular potentials (V e ). One is calculated at a distance 10 μm from the hillock (solid trace) and the other, ∼80 μm from the soma near the apical trunk (dashed trace). The peak-to-peak V e close to the apical trunk is highlighted. D: V e for 30 locations in the x, y plane of the neuron in A. Each trace is 3 ms in duration, and the scale bar is 1 mV. Traces are plotted at the spatial locations at which they were calculated. The soma is centered at (0, 10). E: peak-to-peak amplitude of V e as a function of distance from the hillock. Locations were sampled along the x-, y-, and z-axes at regular intervals. The 1/r best-fit line (least-squares) and equation are shown. Download figureDownload PowerPoint

The polarity of membrane currents differed between neuron sections (Fig. 1B). For example, membrane currents in the soma (and hillock and initial segment) were initially negative, whereas currents in the dendrites were initially positive. As V e is the superposition of potentials generated by all currents, V e could be initially positive or negative depending on location (Fig. 1, C and D). For example, V e calculated close to the soma was dominated by perisomatic currents in the hillock and initial segment (Fig. 1C) compared with V e close to the apical trunk (Fig. 1C). The amplitude of V e had an approximately inverse relationship with distance r, measured from the centroid of the initial segment and hillock (Fig. 1, D and E). For r > 100 μm, peak-to-peak V e was well approximated as 1/r, whereas locations closer than 100 μm had larger-magnitude V e . With the sampling of locations in three dimensions around the initial segment and hillock, the coefficient of best fit of peak-to-peak V e was A e = 7.7 mVμm, such that peak-to-peak V e ≈ A e /r. Thus we note three characteristics of V e : it can be of large amplitude (>1 mV), have both depolarizing and hyperpolarizing phases, and approximate 1/r spatial dependency.

The Spiking Dynamics of Nearby Neurons Can Be Linked by Ephaptic Coupling We next sought to characterize how V e around a spiking neuron affected V m in a nearby neuron. Because V m = V i − V e , ephaptic coupling deflected membrane voltage V m in the receiving neuron with the opposite polarity to V e (Fig. 2, A and B). For this interaction, we did not include ephaptic effects on the stimulating neuron and considered only one-way stimulation. Figure 2A shows V e calculated 5 μm from the stimulating neuron hillock. We sampled the peak-to-peak deflection of V m in the receiving neuron, due to ephaptic stimulation at multiple offsets between the two neurons. Analogous to A e , we computed the coefficient of best fit (A m = 1.5 mVμm), such that peak-to-peak V m ≈ A m /r. Due to the neurons' complex topologies, there was some variation in peak-to-peak V m deflection at a given distance (the deflection in Fig. 2B shows a high-amplitude example for V m deflection at 5 μm), although in general, V m deflection was well fit with a 1/r spatial dependency. Fig. 2.Ephaptic coupling in NEURON. A: stimulating potential V e calculated 5 μm from stimulating neuron hillock. Additional traces were created by multiplying V e by coefficients +0.5 (solid, thin trace), 0 (solid, thick trace), and −0.5 and −1 (dotted traces). B: ephaptic stimulation while receiving neuron is at resting membrane voltage (V m ). Receiving and stimulating neurons were separated by 10 μm. Arrowhead shows start of ephaptic stimulation. C: ephaptic stimulation while receiving neuron is close to threshold. Identical ephaptic stimulation as in B. The arrowhead shows the start of ephaptic stimulation. Up-down V m deflections advance the spike, whereas down-up delay it. Time t = 0 corresponds to the peak of V m with no ephaptic stimulation (nil). D: quantification of C. Data points are generated by multiplying stimulating V e by a scalar coefficient between 1 and −1, as in A–C. E: effect of varying start of ephaptic stimulation (arrowheads in A–C) on spike timing. Both up-down (solid trace) and down-up (dashed trace) V m deflections can advance and delay the spike. Arrowhead corresponds with arrowheads in A–C. Download figureDownload PowerPoint

Single-spike ephaptic coupling altered spike timing in the receiving neuron (Fig. 2C). When a 1.5-mV V e was delivered, ∼1 ms before the peak of the spontaneously generated spike in the receiving neuron, the spike was advanced by close to 0.4 ms. In general, biphasic V e , with an initial negative part occurring at the hillock, led to earlier spike times, whereas an initially positive V e , occurring at the hillock, delayed spike times. The reversal of the polarity of V e at every location of the receiving neuron and the delivery of it at the same time resulted in delayed spike times rather than advanced spike times (Fig. 2D). Due to V e being biphasic, spike-timing effects were also dependent on the time at which V e stimulation was delivered. Changes to V e by an initial negative part could advance spike times, as the initial negative part of V e depolarized the membrane potential and triggered a spike. Biphasic traces of V e with an initial negative and subsequent positive part are here referred to as up-down, as this is the shape of their effect on V m in the receiving neuron. The same up-down stimulation delivered 2 ms earlier had the opposite effect, as the subsequent positive part of V e hyperpolarized V m and delayed the spike time (Fig. 2E). Coupling by a biphasic V e with the opposite polarity (down-up) advanced spike times when delivered a few milliseconds before the spike and delayed spike times when delivered directly before the spike (Fig. 2E).

Global Synchrony and Microclusters Emerge in Small, Ephaptically Coupled Neuronal Networks We extended the NEURON simulations into a network of neurons that were coupled only ephaptically; again, no synaptic or other coupling was implemented. Neurons were organized in a square laminar grid (Fig. 3A). When spacing between neurons was changed, distinct, dynamic phenomena could be observed. At large distances, ephaptic effects were negligible, and neuron firing was asynchronous, as measured by small, pairwise correlations between somatic V m (Fig. 3A). Reduced spacing, e.g., 10 μm, resulted in synchronization between a majority of model neurons (Fig. 3B). Spacing closer than 5 μm resulted in local clusters of synchronous neurons (Fig. 3C). Fig. 3.Ephaptically coupled network in NEURON. A: ephaptically coupled network of 16 neurons arranged in a square, 4 × 4 grid with 5,000 μm spacing. Top: schematic of laminar network. Neurons plotted to scale. Center: gray-scale plot of V m for all 16 neurons. V m color code is shown on the right of the figure; action potential peaks are white lines. Simulations were run in the NEURON environment. Bottom: pairwise correlations among all neurons in the network (Pearson correlation coefficient). B: same as A but for a network with 10 μm spacing. Top: for clarity, only somas are shown (not to scale). C: same as A but for a network with 3 μm spacing. Download figureDownload PowerPoint



Persistent Network Synchrony Requires Homogeneous and Specific Ephaptic Coupling Profiles Because ephaptic potentials are biphasic, an ephaptic potential from one neuron can increase or decrease the phase of another neuron, dependent on their relative phases. We used an analytic model to outline how the shape of the ephaptic potential—up-down or down-up—contributed to synchronization. We first studied two tonically firing neurons (stimulating and receiving neuron, modeled as phase oscillators) with equal intrinsic firing rates. The firing rate of the stimulating neuron was constant and therefore, independent of the activity of the receiving neuron. The receiving neuron received ephaptic potentials from the stimulating neuron, which altered the receiving neuron's phase progression. Specifically, for up-down potentials, a stable steady state existed at a phase difference of zero, whereas for down-up potentials, the stable steady state was nonzero and equal to the zero-crossing phase φ e (see methods for a full description). Thus dependent on the ephaptic potential shape, the two neurons would synchronize or attain a persistent, nonzero-phase difference. For networks of three or more neurons, the shape of the neuron-to-neuron coupling function determines the qualitative phase-difference dynamics on the network level. Specifically, a network with only up-down ephaptic potentials has a network-stable steady state in which all neurons are fully synchronized. In contrast, we mathematically proved that a network with only down-up potentials cannot have a network-stable steady. This finding implies that networks with mostly up-down shape potentials (p+ = 1) would tend toward increased synchrony, whereas networks with mostly down-up potentials (p+ = 0) should undergo accelerated desynchronization. To investigate this dependence of synchrony on p+ further, we determined phase difference distributions and synchrony in networks with N → ∞ neurons and stochastic fluctuations in the phase differences. A network with up-down potentials indeed showed an accumulation close to a phase difference of zero (Fig. 4A). A network with down-up potentials showed a displacement of the distribution away from zero-phase difference, indicating again an active desynchronization (Fig. 4B). In summary, the fundamental shape of the ephaptic potentials introduces an asymmetry in the effects of ephaptic coupling between pairs of neurons on the microscopic level, which feeds into the macroscopic behavior of the neuronal network. For increasing strength of ephaptic coupling, up-down potentials lead to a stronger synchrony increase; down-up potentials increasingly drive desynchronization (Fig. 4C). Fig. 4.General conditions for synchrony. A and B: different effect of up-down vs. down-up ephaptic potentials (c̅ > 0 vs. c̅ < 0, respectively), coupling strength (absolute magnitude c̅), and different temporal profiles (ϕ e ) on synchrony of the phase difference distribution [p(δ)]. c̅ = +5,000 (A); c̅ = −5,000 (B). ϕ e = 0.01, 0.02, … , 0.1 are indicated by shading from black to gray. C: mean z value for different c̅ (ϕ e same as in other panels). Download figureDownload PowerPoint



Emergent Behaviors in Oscillator Networks Depend on Spatial Distribution We assessed the relevance of ephaptic coupling in a spatially distributed network of phase oscillators, which serves as simplified representations of ephaptically coupled neurons. The reduced numerical workload of simulating an oscillator network allowed an extensive assessment of larger networks, which was not feasible in NEURON. Connections between the oscillators mimicked ephaptic coupling in three ways: connections took effect in a temporally distributed manner and were biphasic, coupling amplitude scaled with the distance between oscillators as 1/r, and each oscillator was coupled to every other. Examples of the biphasic ephaptic potentials are shown in Fig. 5A. Their effects on phase progression (Fig. 5C) are defined in Eqs. 11 and 13. Where known, we assigned model parameters in accordance with literature and common values (Table 1). The following parameters were varied between simulations: N, the number of oscillators; R, the distance parameter that sets the density of oscillators and thus also their spacing; and p+, the probability that a given ephaptic connection was up-down (as opposed to down-up). We used the z statistic to indicate population synchrony (z = 0 means full asynchrony; z = 1 means full synchrony). Fig. 5.Network model schematic for 3 coupled neurons. A: geometric arrangement and ephaptic potentials defining the connections. The amplitude of A/r b,c is larger, because r b,c < r a,b . Inset: examples of the biphasic ephaptic potentials are shown. B: oscillator phases without coupling (R → ∞). C: oscillator phases with coupling. Download figureDownload PowerPoint

We found that increasing N and decreasing R led to transitions from asynchronous to synchronous behavior. The transition specifically required a predominance of up-down ephaptic potentials: down-up-dominant and balanced networks did not synchronize (p+ = 0.0, and p+ = 0.5; Fig. 6, A and B), whereas up-down-dominant networks did (p+ = 1; Fig. 6C). This is in agreement with our earlier predictions of how the qualitative shape of the individual ephaptic potentials determines the establishment or loss of network-level synchrony. For intermediate R values, the degree of synchrony was largely dependent on initial conditions (Fig. 6C). Fig. 6.Occurrence of synchrony in the network model. The z statistic of synchrony depends on the number of neurons (N), distance parameter (R), and fraction of up-down potentials (p+). Twelve simulations/data point; z calculated from 2,500 to 3,500 ms; shown are arithmetic means. Initial phases were either uniformly distributed (thin lines) or fully synchronized (thick lines). p+ indicated above A–C. U, uniform distribution. Download figureDownload PowerPoint

Closer investigation of the relevance of initial conditions demonstrated a “critical slowing down” of the network's z dynamics for intermediate R. In systems with critical behavior, typically, two stable states are attainable by subparts of the system. These subparts can successively entrain larger parts of the overall system, thus leading to global dominance of one stable state. Once a stable state has reached global dominance, it can persist for periods significantly exceeding lifetimes of stable states in subparts of the system. In the oscillator network simulation, global dominance of the asynchronous state is indicated by z = 0 and global dominance of the synchronous state by z = 1. For low R, global asynchrony disappears rapidly, and global synchrony persists as a stable global state (Fig. 7A). Subgroups of oscillators could only stably attain synchrony. Accordingly, no critical slowing down occurs. Instead, the entire network rapidly departs from rapid z = 0 initial conditions, frequently attains z = 1 within seconds, or remains at z = 1 initial conditions. For intermediate R, synchronous as well as asynchronous initial conditions persisted for several seconds (Fig. 7B). Synchrony and asynchrony between oscillators are two possible stable states in subparts of the network. Once the entire network was entrained into the synchronous state, emergence of asynchrony became highly unlikely and vice versa. Because the initial conditions were such network-wide, asynchronous (z = 0) or synchronous (z = 1) states, these initial conditions persisted for several seconds. For higher R, synchrony was lost for all initial conditions (Fig. 7C). In this case, synchrony in subgroups of oscillators was not a stable state. Thus only one stable state was possible in such subgroups, and no critical slowing down could be observed. The network rapidly attained z = 0, irrespective of initial condition. Fig. 7.Hysteresis in population synchrony and asynchrony. N = 100 oscillators were started with all ϕ = 0 (gray lines) or with uniformly distributed ϕ (black lines). A: for low R, all simulations develop toward synchrony (z = 1). B: for intermediate R, the initial synchrony level persists. C: for high R, all simulations approach asynchrony (z close to 0). Simulation parameters: p+ = 1; N = 100; square geometry was used. Download figureDownload PowerPoint

Scanning a range of R in more detail showed three major regimes. Regime A: below R ≈ 10 μm, the firing rate was markedly elevated above the oscillators' intrinsic firing rate (≈20 s−1; Fig. 8A). An elevated synchrony (0.5 ≤ z < 1) was present, but full synchrony was not attained. Regime B: for R ≈ 10 μm, the firing rate had dropped to almost the intrinsic firing rate (Fig. 8A), whereas the entire neuronal population was synchronized (Fig. 8B). Regime C: the increase of R beyond 10 μm led to a reduction of z over the course of the simulation. To assess how long a z value taken on at a given point in time persisted after this point, we calculated the autocorrelation time (τ AC ) of the z traces. The computation of the autocorrelation of z and the calculation of the characteristic exponential decay rate of the resulting autocorrelation function yield τ AC , which is a common statistic for analyzing persistence in a time series (Aburn et al. 2012). τ AC increased up to a maximum of ≈15 s, indicating that specific z values persisted for times that exceed the time scale of individual oscillator periods by several orders (Fig. 8C). For R ≥ 80 μm, the average z was close to zero, and τ AC took on a constant value of ∼2.5 s. The dynamics in regimes B and C seem generally in line with our prior reasoning: below a critical R, the population synchronized spontaneously and desynchronized spontaneously above a critical R. An alternative explanation would be that for increasing R, the development of synchrony took increasingly long: below R ≈ 80 μm, synchrony would be attained before the simulation finished; above R ≈ 80 μm, synchrony would have taken longer than the simulation time; for increasing R, the development of global synchrony would have taken longer, leading to an increased τ AC . Because simulating for a longer time could never rule out this interpretation, we executed the same simulations but used perfectly synchronous initial conditions. We found that regime A also exhibited an increased firing rate and intermediate z values, as observed before (Fig. 8, D and E). In regime B, synchrony was maintained throughout the simulation, leading to z = 1, whereas synchrony was lost spontaneously throughout regime C (Fig. 8E). The loss of synchrony was more pronounced for increasing R, indicating a more rapid decrease of z after the simulation is initialized. Lastly, a slight increase in τ AC was visible when the transition from regime B to regime C was approached (Fig. 8F). This indicates occasional reductions of z throughout the simulation, which occurred on a somewhat slower time scale than the individual neurons' intrinsic firing rate. Fig. 8.Influence of distance parameter R on global synchrony. N = 200 oscillators with only up-down potentials (p+ = 1) was initialized with random phases (A–C) or synchronized phases (D–F). Statistics calculated are average frequency of individual neuron spikes (f; A and D; gray lines indicate intrinsic frequency of 20 s−1), synchrony (z; B and E), and the autocorrelation time of the z time course (τ AC ; C and F). An elongated geometry was used; simulations run for up to N·6,000 iterations or until 40 s elapsed time; statistics calculated for t ≥ 2 s. Gray boxes, individual simulation runs; solid lines, arithmetic mean of individual simulations. Download figureDownload PowerPoint

In summary, we found three R-dependent regimes: (A, subcritical regime) markedly accelerated firing rate and intermediate synchrony; (B, critical regime) slightly accelerated firing rate, spontaneous emergence of synchrony across the entire population, and persistence of states of synchrony or asynchrony over extended periods of time; and (C, supercritical regime) firing at individual neurons' intrinsic firing rate and spontaneous loss of synchrony. Regimes B and C are explicable by our prior reasoning. A network-stable, synchronous state was established and lost dynamically, which became more pronounced for larger amplitudes of the ephaptic potentials. In regime A, however, we noticed indications of short-range synchronization, which could not be assessed using a global measure of synchrony, such as the z value.