## Abstract

Neuronal hypersynchrony is implicated in epilepsy and other diseases. The low-frequency, spatially averaged electric fields from many thousands of neurons have been shown to promote synchrony. It remains unclear whether highly transient, spatially localized electric fields from single action potentials (ephaptic coupling) significantly affect spike timing of neighboring cells and in consequence, population synchrony. In this study, we simulated the extracellular potentials and the resulting coupling between neurons in the NEURON environment and generalized their connection rules to create an oscillator network model of a sheet of ephaptically coupled neurons. With the use of both models, we explained several aspects of epileptiform behavior not previously modeled by synaptically coupled networks. Importantly, reduction of neuron spacing induced synchronization via single-spike ephaptic coupling, agreeing with seizure suppression seen clinically and in vitro via extracellular volume adjustment. Further reduction of neuron spacing yielded locally synchronized clusters, providing a mechanism for recent in vitro observations of localized neuronal synchrony in the absence of synaptic and gap-junction coupling.

- ephaptic coupling
- synchrony
- mathematical model
- oscillator network

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.

## METHODS

### NEURON Simulations

#### Ephaptic coupling between two neurons.

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.

At location *x*, the *V*_{e} (*x*, *t*) was calculated as the sum of individual potentials (*v*_{i}) generated by all membrane currents
(1)
where *v*_{soma} is the potential due to somatic currents (*I*_{soma}). Since the soma is assumed to be spherical, *v*_{soma} can be calculated as
(2)
where *d* is the distance between *x* and 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 (*v*_{i}) were calculated for each cylinder. For full details, see Holt (2008).

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).

#### Ephaptically coupled, small networks in NEURON.

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}.

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.

### Mathematical Analysis of Conditions for Global Synchrony

#### Phase relationship between two neurons.

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
(3)
where *f* is *δ*-periodic on the interval [0, 1). *f* was specified as a function of *δ* based on qualitative features of our NEURON results to allow general mathematical conclusions. Quantitative parameters were assumed only in later simulations of network dynamics.

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
(4)

#### Conditions for a network-stable synchronous state.

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} = = 0; all phase differences vanish, and the network-stable state is perfectly synchronous. For the opposite case, *p*^{+} = 0, all *c* < 0, so
(5)
for all *a* ≠ *b* is the condition for a global steady state. For any two oscillators (*a* and *b*) from the network,
(6)
must hold to fulfill the condition from *Eq. 5*. With the introduction of a third oscillator (*c*), the possibility of a network-wide steady state can be disproved by the following contradiction with *Eq. 6*
(7)

#### Distribution of phase differences in an ensemble of oscillators.

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)
(8)
where *D*, the diffusion constant, quantifies the strength of fluctuations in *δ*. At steady state
(9)
where *c̄* = *c*/*D*. The parameters *c* and *D* are thus described by the single parameter *c̄*. By separation of variables
(10)

Σ is a normalization constant that ensures ∫_{0}^{1} *p**(*δ*)*dδ* = 1. To compute *p**(*δ*), an approximated piecewise linear *f*(*δ*) function is used (see Fig. 5*A*)
(11)
The *z* statistic, commonly used to quantify synchrony of phase oscillators, then is(12)

### Ephaptic-Like Oscillator Network Model

*N* phase oscillators representing ephaptically coupled neurons were placed at random positions in a square with side-length *R*. *R* effectively defines the oscillator density, *N*/*Area* = *N*/( *R*)^{2} = 1/*R*^{2}. The phase *ϕ*_{i} ∈ [0, 1) of neuron *i* progresses according to
(13)
with *i*, *j* = 1, . . . , *N*. The intrinsic frequencies (*ω*_{i}) were randomly drawn from a normal distribution centered on *ω*_{0} with an SD of *σ*_{ω}. The other terms represent ephaptic coupling of other neurons *j* to neuron *i*, which were summed to give the overall alteration of *ϕ̇*_{i} by ephaptic coupling. *r*_{ij} is the Euclidean distance between neurons *i* and *j*. The polarity of the ephaptic-like coupling was set by *ĉ*_{j,i} = ±1 (i.e., *ĉ*_{j,i} = +1 for “up-down” or *ĉ*_{j,i} = −1 for “down-up”; *ĉ*_{j,i} = 0 was used when *j* = *i* to prevent oscillators from coupling to themselves). *ĉ* is analogous to parameters *c* and *c*̄ in that all determine the polarity of ephaptic coupling, although *ĉ* does not control the amplitude of coupling. *ĉ*_{j,i} was assigned randomly for each pairwise connection, with a probability *p*^{+} that *ĉ*_{j,i} = +1 (e.g., *p*^{+} = 1 means that all ephaptic potentials have an up-down shape). *H* (*ϕ*_{i}/*ω*_{i} − *T*_{r}) is a Heaviside function, ensuring that the oscillator *i* was not affected by ephaptic potentials during a refractory period of duration *T*_{r} after crossing *ϕ*_{i} = 1. Outgoing ephaptic potentials were triggered when *ϕ*_{i} = 1 was crossed. At the same time, *ϕ*_{i} was reset to zero. The three parameters *N*, *R*, and *p*^{+} were changed among our different simulations. The fixed simulation parameters are shown in Table 1.

The detailed effect of ephaptic potentials of oscillator *j* on oscillator *i* was described by
(14)
*T*_{e} is the time at which the ephaptic potential deflection first returns to 0. The coefficient *A* expresses 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 the *layer V* pyramidal neuron simulated in NEURON, this corresponds to an ≈20-mV difference in membrane potential (Table 1). Therefore, the amplitude of the change in *ϕ* induced by an ephaptic potential at 1 μm was *A*_{m}/20 mV. *A*_{m} describes the peak-to-peak deflection of *V*_{m} when 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, *T*_{e}/2. Thus
(15)

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.

Population synchrony was calculated with the *z* statistic,
16
where *z* = 1 means complete synchrony (all phases *ϕ*_{j} are equal), and *z* = 0 means complete asynchrony (uniform distribution of phases).

### Simulation of Experiments

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*).

Hierarchical clustering (single-linkage algorithm) of the oscillators was executed based on the average minimum distance (AMD) between spike trains (spikes detected at *ϕ* = 1 threshold) (Feldt et al. 2009). The AMD between two oscillators *m* and *n* was calculated as
(17)
where Δ*t*_{k}^{m,n} is the smallest time difference from the time of spike *k* in spike train *m* to any spike in spike train *n*, and *N*_{m} is the number of spikes detected in spike train *m*. The number of clusters was chosen by maximal modularity using the similarity measure
(18)

## 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. 1*A*). 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.

The polarity of membrane currents differed between neuron sections (Fig. 1*B*). 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. 1*C*) compared with *V*_{e} close to the apical trunk (Fig. 1*C*).

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 2*A* 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. 2*B* 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.

Single-spike ephaptic coupling altered spike timing in the receiving neuron (Fig. 2*C*). 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. 2*D*).

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. 2*E*). 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. 2*E*).

### 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. 3*A*). 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. 3*A*). Reduced spacing, e.g., 10 μm, resulted in synchronization between a majority of model neurons (Fig. 3*B*). Spacing closer than 5 μm resulted in local clusters of synchronous neurons (Fig. 3*C*).

### 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. 4*A*). A network with down-up potentials showed a displacement of the distribution away from zero-phase difference, indicating again an active desynchronization (Fig. 4*B*). 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. 4*C*).

### 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. 5*A*. Their effects on phase progression (Fig. 5*C*) 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).

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. 6*C*). 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. 6*C*).

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. 7*A*). 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. 7*B*). 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. 7*C*). 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.

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. 8*A*). 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. 8*A*), whereas the entire neuronal population was synchronized (Fig. 8*B*). *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. 8*C*). 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. 8*E*). 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. 8*F*). This indicates occasional reductions of *z* throughout the simulation, which occurred on a somewhat slower time scale than the individual neurons' intrinsic firing rate.

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.

### Global and Microcluster Synchrony Emerge in an Elongated Neuronal Field When Neuron Density Is Increased

We adapted our oscillator model to the elongated geometry of neuronal fields seen in recent experiments with mouse hippocampal slices. In this geometry, we assessed the development of synchrony in microclusters for different values of the distance parameter *R* (Muldoon et al. 2013). For *R* values below those that show *z* ≈ 1 (subcritical), the neuronal field separated into localized clusters of synchrony; transient synchrony between the local clusters occurred (Fig. 9, *A* and *B*). For *R* values that show *z* ≈ 1 (critical), global synchrony throughout the neuronal field was seen (Fig. 9*C*); clusters spanned large parts of the neuronal field and were not as clearly separated any longer (Fig. 9*D*). For *R* values exceeding the range in which *z* ≈ 1 (supercritical), no clear pattern of global or local synchrony emerged (Fig. 9*E*); correspondingly, random cluster associations occurred, which spanned the entire neuronal field with no apparent spatial order (Fig. 9*F*). A more systematic assessment for different values of *R* confirmed these observations as general patterns (Fig. 9*G*). For low *R* values (*R* < 10 μm), few clusters (<10) occurred, which spanned only small fractions of the neuronal field's length. For increasing *R* (10 μm < *R* < 40 μm), even fewer clusters occurred (<5), which spanned a large part of the neuronal field's length, as was the case when global synchrony was established within the neuronal field. For higher *R* values (*R* > 40 μm), the number of clusters increased, and the clusters spanned large parts of the neuronal field, as observed in Fig. 9*F*. Exclusively for *R* values that led to microclustering, the average firing rate within a cluster increased with the number of oscillators within a given cluster (Fig. 9*H*).

## DISCUSSION

### Functional Relevance of Single-Spike Ephaptic Coupling

We demonstrated that single-spike ephaptic coupling could significantly affect spike timing between pairs of model neurons, both in the relatively realistic NEURON environment and a simplified oscillator network. *V*_{e} around a spiking neuron should, by the laws of electromagnetism, affect charged ions and the activity of nearby neurons. However, the effects are generally low magnitude (*V*_{e} < 1 mV), transient (<10 ms), and due to the radial symmetry of *V*_{e}, nonspecific. For these reasons, the functional relevance of single-spike ephaptic coupling is questioned (Buzsáki et al. 2012). Despite such reservations, recent experiments performed by Anastassiou et al. (2011) and Frohlich and McCormick (2010) showed that sub-millivolt changes to *V*_{e} can increase synchronization. These stimulating potentials were of similar or lower magnitude to those produced around a single spiking pyramidal neuron (Gold et al. 2006; Holt and Koch 1999).

Aside from Park et al. (2005), most previous studies of cortical ephaptic coupling used low-frequency (e.g., 1 s^{−1}), LFP-like stimulating potentials (Anastassiou et al. 2011; Berzhanskaya et al. 2013; Frohlich and McCormick 2010; Radman et al. 2007). However, neurons are also subject to higher frequency extracellular voltages from single spikes. Intracellular recordings of mouse hippocampal slices show that the fields around a single-spiking neuron cause “spikelets,” biphasic voltage deflections, in intracellular recordings of nearby neurons that last, at most, 10 ms (Vigmond et al. 1997). Whereas these spikelets are indeed transient, we showed here that *V*_{e} of their magnitude and duration could nevertheless affect spike timing when the receiving neuron is close to threshold (Fig. 2) and synchronize ephaptically coupled model networks (Figs. 3 and 6–9). For symmetric, biphasic spikelets, as in our model, the effect of ephaptic coupling is thus dependent on the phase of the receiving neuron being close to threshold. This phase dependence increases the effective specificity of single-spike ephaptic coupling. Areas, such as the hippocampus, have relatively high neuronal densities, increasing the likelihood that a neuron in its “receptive” phase (close to threshold) is in the relevant spatial proximity.

### Relevance of Profile and Distribution of Ephaptic-Like Potentials

Our simplified oscillator network model and simulated networks within NEURON demonstrated that single-spike ephaptic coupling could lead to increased synchronization. Additionally, we made testable predictions about the role of the shape of *V*_{e} (up-down vs. down-up) and how synchronization and neuronal clustering vary as a function of neuronal density. Since a predominance of up-down-type ephaptic coupling (*p*^{+} ≈ 1) was required for the occurrence of synchrony (Fig. 6), we wondered if different geometric arrangements of neurons are more likely to yield *p*^{+} ≈ 1 and thereby promote synchrony. In our model of a *layer V* pyramidal cell, the somatic region (soma, hillock, initial segment) produced conditions for *p*^{+} = 1 via large contributions from Na^{+} currents (initially, negative *V*_{e}), followed by slower K^{+} currents (subsequently, positive *V*_{e}). Fields around the dendrites had opposite polarity (*p*^{+} = 0). However, since *V*_{e} around the soma was at least an order of magnitude greater than dendritic *V*_{e}, significant ephaptic coupling strengths were limited to a zone of *p*^{+} = 1 around the soma.

Although the shape of the extracellular potential is dependent on the site of action-potential initiation (Linden et al. 2010), empirical recordings and simulations suggest that the perisomatic region is the dominant location for *p*^{+} ≈ 1 (Buzsáki et al. 2012; Gold et al. 2006). A laminar organization, as seen in the cortex and hippocampus, aligns neurons' *p*^{+} = 1 zones and excludes *p*^{+} < 1 factors from this region, thereby promoting ephaptic synchronization. Other properties of pyramidal cells, such as their elongate shape that can span up to 1 mm, make them more sensitive to extracellular potentials compared with other cell types (Radman et al. 2009). Therefore, pyramidal cells organized in a laminar structure, such as in the hippocampus, may be especially susceptible to synchronization from ephaptic coupling.

### Ephaptic-Like Coupling Can Explain the Emergence of Synchronized Microclusters

For the lowest *R* values that we investigated in our oscillator network, ephaptic-like coupling was strong enough to allow the establishment of locally synchronized microclusters (Figs. 3 and 9). Within microclusters, the firing frequency was higher in bigger clusters. Also, clusters transiently aligned their firing patterns with those of other clusters. In dentate gyrus slices from epilepsy model mice, where synaptic and gap-junction coupling was blocked, the same observations were made (Muldoon et al. 2013). Furthermore, the neurons were driven into spontaneous firing by a high extracellular potassium concentration ([K^{+}] = 7.5 mM), matching our model assumptions. We also saw locally synchronized clusters of ephaptically coupled neurons within our NEURON network simulations. These clusters arose spontaneously from asynchronous initial conditions. Whereas the mechanism of microcluster formation was unclear in the experimental study of Muldoon et al. (2013), we demonstrated here that ephaptic coupling is sufficient to explain many details of the experimental observations.

### Neuronal Spacing and Epileptogenesis

Dudek et al. (1990) and Roper et al. (1992) suppressed all synaptic communication in rat hippocampal slices consisting of pyramidal neurons, after which, they proceeded to control the extracellular volume via control of extracellular osmolality by membrane-impermeant solutes or low osmolality medium. Reduction of extracellular volume (corresponding to *R* reduction in our oscillator model) reproducibly increased the frequency of occurrence of epileptiform dynamics, whereas an increase of the extracellular volume (*R* increase) abolished the occurrence of epileptiform dynamics. These findings establish that under conditions of suppressed synaptic coupling in laboratory experiments, spacing between neurons (as altered by volume compression and expansion) controls the occurrence of epileptiform dynamics on the level of neuronal networks. Our model results, based solely on ephaptic coupling, demonstrate a sufficient and specific mechanism for these observations.

The findings of Dudek et al. (1990) and Roper et al. (1992) do not answer whether synaptic volume compression-independent mechanisms take precedence over these nonsynaptic mechanisms in humans under live conditions and how this impacts actual epileptic episodes. Haglund and Hochman (2005) investigated a group of patients suffering from cortical and mesical temporal-lobe epilepsy for their reaction to furosemide and mannitol administration, which increase water content and volume of the extracellular fluid (corresponding to *R* increase in our model). In agreement with the earlier laboratory findings, mannitol and furosemide administration suppressed the spontaneous occurrence of electrical spikes, the capability of external induction of epileptiform discharges, and the spread of epileptiform activity.

### Hysteresis of Ephaptically Mediated Synchrony

Our simulations showed a critical *R* regime where the asynchronous and the synchronous states persist for prolonged periods once they are established. This could explain experimental observations of prolonged epileptiform and desynchronized phases without the need to invoke additional, slower processes, such as changes in the extracellular fluid composition (Jiruska et al. 2010a).

### Model Limitations

We note a number of limitations in our NEURON and oscillator treatments. First, our NEURON simulations were based on a single type of neuron. It is therefore not clear how far our findings can be generalized to other neuron types. Given that ephaptic coupling mostly advances or delays existing spikes, it is likely that results here are dependent on similar intrinsic firing rates among model neurons. Second, we assessed synchronization via ephaptic coupling and disregarded all other modes of coupling. This reflects the conditions of the experiments explained by our model. This does not, however, address the dynamic interplay between ephaptic and other modes of coupling (e.g., synapses and gap junctions) in in vivo neuronal networks. Third, the small population of neurons and the simplified two-dimensional rectangular architecture fall short of the geometric complexity of real neuronal architectures. Finally, our model assesses ephaptic coupling only as a physical mechanism, not as an aspect of an information-processing system. It remains largely unclear how ephaptic coupling and the resulting synchrony integrate with the information processing by neuronal networks.

We used extracellular potentials with a large amplitude—our value for *A*_{e} may be larger than in other simulations [e.g., Gold et al. (2006)]. The *layer V* pyramidal cell that we model likely produces extracellular potentials with higher than average amplitude. It is also possible that experimental measures of extracellular voltage are artificially lowered by scarring around the electrode. Furthermore, both the positive and the negative part of *f* have the same magnitude, whereas this is not the case for most extracellular action potentials (Gold et al. 2006).

### Comparison with Existent Models

The piecewise linear approximation used in the oscillator network model captures the biphasic shape of the extracellular potential but not the detailed nuances of the potentials and their dependence on relative neuron position. Together with the decision to model neurons as simple phase oscillators, this is a significant reduction of complexity and realism compared with the physically realistic NEURON simulations. This reduction in complexity made it feasible to treat larger networks. In spite of the reduced realism, our phase-oscillator models demonstrated several principles of ephaptic coupling as a mechanism for neuronal synchrony and explained several experimental observations in sufficient mechanistic detail. The observation of global synchrony and synchronized microclusters at different densities of oscillators was also confirmed in NEURON simulations.

A previous study by Park et al. (2005) provides complementary insights, addressing the coupling functions between a pair of neurons in more detail. Park et al. (2005) first simulated a pair of two-compartment neuron models that are coupled via ephaptic and synaptic connections to derive detailed coupling functions for a phase-oscillator model of the pair of neurons. Accordingly, their approximation successfully explains detailed phase-locking patterns observed in the compartment model (Park et al. 2005). Thus for the long-term dynamics of two coupled neurons, Park et al. (2005) give a more detailed account.

A number of mathematical models of epileptiform behavior are based on changes in synaptic connections [see Lytton (2008) for review]. In one class of models, “rewiring” of the normal synaptic network is assumed as the main mechanism. Positive-feedback connections (potentiating neuronal activity) are added to the network, inhibitory connections in the network (suppressing excessive activity) are removed, or the neuronal network topology is altered. For example, Dyhrfjeld-Johnsen et al. (2007) constructed a realistic “virtual dentate gyrus” model. Increasing sclerosis of the dentate gyrus was mimicked by *1*) increasing removal of neurons with long-range connections and *2*) increasing “sprouting,” the introduction of new, synaptic connections. With increasing sclerosis, network topology and dynamics underwent a biphasic change. Intermediate degrees of sclerosis increased long-range connectivity (“small world” network topology) and led to hyperexcitability of the network, thus supporting epileptogenesis. Severe sclerosis, however, removed long-range connections to such an extent that it broke down effective long-range connectivity, thus preventing long-range spread of activity and hyperexcitability. In another realistic dentate gyrus model, based on synaptic connections that were investigated by Morgan and Soltesz (2008), positive-feedback connections, specifically within hubs of granule cells, established epileptogenic microcircuits. Similarly, a simplified model, based on a small network constructed by Netoff et al. (2004), could change from quiescence to excitability to epileptogenic behavior by network structural changes. Both models by Morgan and Soltesz (2008) and Netoff et al. (2004) predicted traveling waves of neuronal hyperactivity. In the experiment, however, spatially stationary clusters of synchrony and hyperactivity are observed (Jiruska et al. 2010a; Muldoon et al. 2013).

Berzhanskaya et al. (2013) constructed small networks representing hippocampal layers using NEURON. Significant geometric realism was considered in constructing these networks, and pyramidal neurons as well as interneurons were included in the simulation. Model assumptions regarding specific ionic currents, synapse distribution, and the response to external electrical fields were assessed in single neuron experiments (Berzhanskaya et al. 2013). This work therefore addresses many of the simplifications that we made in our study. However, ephaptic coupling was not considered.

Importantly, experimental studies show that with all synaptic connections silenced, seizure-like synchrony can exist in hippocampal slices (Jefferys and Haas 1982). Thus other mechanisms, such as ephaptic coupling, are likely involved in epileptiform behavior (Jefferys 1995; Weiss and Faber 2010). To our knowledge, our study represents the first network model investigating synchronization solely via single-spike ephaptic coupling.

### Future Directions

Additional NEURON simulations realistically implementing extracellular vs. intracellular volume changes would deepen our understanding of the connection between volume changes and epileptiform dynamics in neuronal populations. We interpreted *R*, a primary parameter governing the existence of global and cluster synchrony in our simulations, as an effective distance between close-by neurons' membranes. Presumably, changes in the extracellular vs. intracellular volumes would alter the distance between neurons' membranes, which in turn, alters *R*. Following this assumption, we could explain experimental and clinical observations. However, future simulations should be executed where the geometry of the cell membrane is altered to assess this hypothesis fully.

In vitro experiments with slices of neuronal tissue that showed global and cluster synchrony can be extended to assess the effect of extracellular vs. intracellular volume changes (Muldoon et al. 2013). Global and cluster synchrony in these experiments was governed by two parameters: addition of positive ions increased the general neuronal activity, whereas an epilepsy-disposing mutation induced cluster synchrony (Muldoon et al. 2013). Our simulations, along with previous experimental and clinical results (Dudek et al. 1990; Haglund and Hochman 2005; Roper et al. 1992), suggest that the extracellular vs. intracellular volume (adjustable via, e.g., sugars that increase medium osmolality) would exert additional control over global and clustered synchrony.

Lastly, our study suggests that synchrony caused solely by ephaptic coupling will be less likely where neurons are organized less rigidly into a laminar sheet. Synchrony is often studied in tissues with laminar organization of neurons, because only here can EEGs detect synchrony effectively. However, the imaging of individual neurons in a larger field could confirm or falsify the hypothesis that more irregularly structured neuronal tissues are less likely to develop synchrony solely due to ephaptic coupling.

## GRANTS

Funding was for this study was provided by Natural Sciences and Engineering Research Council of Canada (NSERC); McGill Faculty of Medicine Scholarship (to R. G. Stacey); Canadian Institutes for Health Research (CIHR) McGill Systems Biology Graduate Training Program (to L. Hilbert); and Heart and Stroke Foundation (to T. Quail).

## DISCLOSURES

The authors declare no competing financial interests.

## AUTHOR CONTRIBUTIONS

Author contributions: R.G.S., L.H., and T.Q. conception and design of research; R.G.S., L.H., and T.Q. analyzed data; R.G.S., L.H., and T.Q. interpreted results of experiments; R.G.S. and L.H. prepared figures; R.G.S. and L.H. drafted manuscript; R.G.S., L.H., and T.Q. edited and revised manuscript; R.G.S. and L.H. approved final version of manuscript.

## ACKNOWLEDGMENTS

The authors thank Costas Anastassiou for an excellent literature overview, Daniel Zysman for comments on our manuscript, and Michael Mackey for his guidance in general and insisting on physical realism specifically.

## Glossary

- δ
_{a,b} *Difference in phase between oscillator*a*and oscillator*b- Δt
_{AP} *Ephaptically driven change in spike time, measured at peak*V_{m},*s*- ρ
_{e} *Bulk extracellular resistivity, Ωm*- σ
_{ω} *SD of oscillator frequencies in an oscillator network, s*^{−1}- ϕ
_{e} *Zero-crossing phase of*f,*ms*- ϕ
_{i} *Phase of oscillator*i- ω
_{i} *Frequency of oscillator*i,*s*^{−1}- ω
_{0} *Mean frequency of all oscillators in an oscillator network, s*^{−1}- A
*Amplitude of ephaptic-like coupling, μm s*^{−1}- A
_{e} *Peak-to-peak extracellular voltage*V_{e}*at 1 μm, mVμm*- A
_{m} *Peak-to-peak deflection of membrane voltage*V_{m}*at 1 μm, mVμm*- c
*Controls the amplitude and polarity of ephaptic-like potential*- ĉ
*Controls the polarity of ephaptic-like potential*- c̄
*Parameter combining*c*and*D- d
*Distance between*x*and soma center for NEURON calculations, μm*- D
*Diffusion constant quantifying strength of fluctuations in*δ- d
_{m,n} *Average minimum distance, used to cluster spike trains*- f
*Ephaptic-like potential*- g
*Function describing the sign of ephaptic-like coupling, s*^{−1}- I
_{soma} *Membrane currents in the soma, A*- N
*Number of oscillators in an oscillator network*- p
^{+} *Probability of up-down ephaptic-like potential*- p(δ)
*Probability distribution of phase differences*- R
*Parameter controlling spacing (density) of oscillators in oscillator network, μm*- r
_{ij} *Distance between oscillator*i*and oscillator*j,*μm*- T
_{e} *Zero-crossing time of*g,*s*- T
_{r} *Refractory period of oscillators, s*- V
_{e} *Extracellular potential, mV*- v
_{i} *Contribution to*V_{e}*from straight line segment*i,*mV*- V
_{m} *Membrane potential, mV*- v
_{soma} *Contribution to*V_{e}*from soma, mV*- x
*Location in NEURON model, μm*- z
*Synchrony statistic*

- Copyright © 2015 the American Physiological Society