A Role for Electrotonic Coupling Between Cortical Pyramidal Cells

Many brain regions communicate information through synchronized network activity. Electrical coupling among the dendrites of interneurons in the cortex has been implicated in forming and sustaining such activity in the cortex. Evidence for the existence of electrical coupling among cortical pyramidal cells, however, has been largely absent. A recent experimental study measured properties of electrical connections between pyramidal cells in the cortex deemed “electrotonic couplings.” These junctions were seen to occur pair-wise, sparsely, and often coexist with electrically-coupled interneurons. Here, we construct a network model to investigate possible roles for these rare, electrotonically-coupled pyramidal-cell pairs. Through simulations, we show that electrical coupling among pyramidal-cell pairs significantly enhances coincidence-detection capabilities and increases network spike-timing precision. Further, a network containing multiple pairs exhibits large variability in its firing pattern, possessing a rich coding structure.


INTRODUCTION
Synchronized neuronal activity in the cortex is essential for information processing underlying many cognitive functions, such as learning, attention, and memory formation (Wang and Buzsaki, 1996;Buzsaki and Draguhn, 2004;Wang et al., 2010). Several sensory regions in the brain, including the olfactory, visual, and auditory cortices, were observed to exhibit synchronized spiking activity when presented with a stimulus (Singer and Gray, 1995;Moore et al., 2001;Haddad et al., 2013).
Experimentalists and computational neuroscientists have shown that, in addition to neuronal communication via chemical synapses, electrical coupling among interneurons, local inhibitory neurons, plays an essential role in generating and maintaining synchronous activity among neurons in the cortex (Gibson et al., 1999;Beierlein et al., 2000;Chow and Kopell, 2000;Tamas et al., 2000;Traub et al., 2001;Amitai et al., 2002;Nomura et al., 2003;Bennett and Zukin, 2004;Ostojic et al., 2009). Such coupling occurs through protein channels called gap junctions that directly connect the interior contents of the cells (Revel and Karnovsky, 1967).
Yet, cortical networks have been shown to exhibit synchronized oscillations typically attributed to gap-junction connectivity even when the gap junctions among interneurons are blocked, suggesting a possible existence of electrical coupling among the excitatory, pyramidal cells (Mercer et al., 2006;Dere and Zlomuzica, 2012). Gap-junction connections between pyramidal cells in the olfactory system of Drosophilia have been shown not only to exist but to be crucial for odor-evoked lateral excitation, a mechanism by which the organism responds to synchronized, odor-specific input (Wilson, 2013). Experiments have shown evidence supporting the existence of electrical connections between the axons of pyramidal cells in the hippocampus (Schmitz et al., 2001), while computational studies have shown that this network of axo-axonal electrically-coupled pyramidal cells can exhibit fast-frequency oscillations and sharp wave-ripple oscillations (Traub et al., 1999(Traub et al., , 2012. Another computational study showed that vast networks of electrical coupling among pyramidal cells can produce spatiotemporal patterns of activity in the cortex (Traub et al., 2010).
There has been little experimental evidence, however, that gap-junction coupling exists between pyramidal cells in the mammalian neocortex. Mercer et al. (2006), while measuring gap-junction connectivity in the rat neocortex, discovered one connected pair of pyramidal cells, and later, in 2010, Wang et al. (2010) measured ten pairs of pyramidal cells in the prefrontal and visual cortices of rats and ferrets to have pair-wise electrical coupling deemed "electrotonic" coupling to distinguish it from the typical gap-junction coupling among interneurons.
The electrotonic couplings (ECs) were measured to be rare, occurring only pair-wise between cells that have touching or overlapping soma, with a coupling probability of 5% (Wang et al., 2010), significantly lower than those assumed in previous computational studies and much lower than the 60% coupling probability observed among interneurons separated up to 80 microns Hestrin, 1999, 2001;Gibson et al., 1999). The junctional conductance of ECs was measured as nearly 25times higher than that of gap junctions between interneurons. Moreover, action potentials were shown to propagate through the junction, a property that has never been measured for gapjunction-coupled interneurons. The protein that might form this strong, but rare, electrotonic coupling, however, remains unknown, making further experimental investigation difficult and network effects of such a junction undetermined (Connors and Long, 2004;Mercer et al., 2006).
In this work, we address a possible role that EC between pairs of pyramidal cells could play in shaping the dynamics of neocortical networks. Through constructing and simulating a realistic model network, we propose that pairs of electrotonicallycoupled pyramidal cells, or electrotonic pairs, act as sensitive coincidence detectors, eliciting synchronized activity in the network in response to coincident input. Further, we describe how electrotonic pairs could evoke many synchronized network events, with each event more tightly synchronized than those elicited by pyramidal cells without EC. Our results further suggest that multiple electrotonic pairs evoke network events which interact with one another, resulting in a network that can exhibit high variability in its spiking patterns, and thus possesses a rich coding structure.

METHODS
To elucidate the biological function of rare electrotonic pairs in a downstream patch of cortical cells, we set up a multi-layer model network in which the electrotonic pair(s) receive varying levels of synchronized input from an upstream network. Neurons in the downstream network are modeled using the Hodgkin-Huxley (HH) neuronal equations and are organized on a grid, with biologically-relevant connection probabilities and strengths for both synaptic and electric connections and for populations of both inhibitory (interneurons) and excitatory (pyramidal) cells, see below for details. Each neuron, including the electrotonic pair, receives a Poisson background drive modeling incoming spikes from neighboring regions. Since realistic features of network dynamics during sensory stimuli may not be adequately captured by Poisson statistics, we simulate an upstream network using conductance-based Integrate-and-Fire (IAF) neurons. The output of this IAF network is used to mimic sensory information transmitted to the electrotonic pair(s) in the downstream network. The activity and synchrony of the downstream patch are analyzed as a function of the input synchrony, or amount of coincident spikes, received by the electrotonic pair(s).
The following subsections discuss the details of the mathematical models of both the downstream model network and upstream input network. A description of the methods used to quantify synchrony is included in this section as well.

Mathematical Models
The IAF Upstream Model The dynamics of the ith model neuron in the IAF upstream network satisfy the set of equations where v i (t) is the voltage, or membrane potential, of the ith neuron, g L is the leak conductance, and g Q i (t) is the synaptic conductance represented by the characteristic shape of experimentally-observed conductance traces (Johnston and Brown, 1983). The right-hand side of the dynamical equation for h Q i contains two sums: The first sum models incoming spikes from the jth presynaptic neuron at times T k j with strength S Q j for both excitatory, Q = E, and inhibitory, Q = I, synaptic inputs. The second sum models incoming spikes to neuron i at times T k i with strength f Q that originate from outside the model network. These external spikes are modeled by a Poisson spike train with rate ν. The dynamics of the voltage are such that the incoming spikes from other neurons in the network and inputs from the Poisson spike train increase (excitatory input) or decrease (inhibitory input) the voltage and, in the absence of incoming spikes, the voltage decays exponentially toward the resting potential, ǫ R . If the voltage is raised such that it reaches a threshold, determined by ǫ T , the neuron is said to have spiked, the spike time is recorded, and the voltage is set to the reset value, ǫ R . The procedure for calculating the spike time and efficiently implementing the time-evolution equations can be found in Shelley and Tao (2001). The parameter values used in this model are the dimensional versions of those used in Shelley and Tao (2001) and can be found in Table 1. The IAF model network contains 100 neurons that are allto-all connected, including 75% excitatory and 25% inhibitory neurons, and results were obtained from simulations of 5 s for each trial.

The HH Downstream Model
The neurons in the downstream network are modeled using a modified version of the HH equations, with an additional ohmic term describing the current through the junctions, for both interneurons and pyramidal-cell pairs. Previous work has modeled gap junctions among interneurons using IAF neurons with the inclusion of an instantaneous jump in the voltage of the post-junctional cell in response to spiking in the pre-junctional cell (Lewis and Rinzel, 2003) or by inserting characteristic action potentials using a spiking kernel (Chow and Kopell, 2000). Here, we use the more complex HH model because EC has a very high conductance, strongly coupling the membrane potential of the pair of neurons at all times, including during an action potential. As was shown experimentally by Wang et al. (2010), EC can propagate full action potentials, producing changes in the voltage of the post-junctional cell that depend on the shape and size of the pre-junctional action potential. The voltage of the ith downstream neuron is described by the equation whereḡ Na andḡ K are the maximal sodium and potassium conductances, and v Na and v K are the sodium and potassium reversal potentials, respectively. In addition, activation and inactivation parameters, m, h, and n, modeling the opening and closing of the ion channels, exhibit dynamics described by the differential equation with each rate variable described by a set of voltagedependent equations as determined by Pospischil et al. (2008) using fitting techniques to match both fast-spiking (interneurons) and regular-spiking (pyramidal) cells' activity in the cortex.
The synaptic conductance is described using fourth-order kinetics as was done by Sun et al. (2010) by the equations Note that by using a smooth function h(v), the synaptic interactions are no longer event-driven, but instead depend on the value of the voltage of the presynaptic neuron. The second term in Equation (1) sums over all pre-synaptic neurons j, which are chosen randomly with probabilities as described in Table 2. The parameter set used in this work, chosen
through matching with voltage-clamp experimental data for electrotonically-coupled pyramidal cells (Wang et al., 2010) and gap-junction connected fast-spiking neurons (Galarreta and Hestrin, 1999), can be found in Table 2. The networks consists of 75% excitatory (pyramidal) cells and 25% inhibitory cells as is understood to occur in the neocortex of mammals (Beaulieu, 1993) and has been done in previous realistic modeling studies (Cai et al., 2005).

Network Organization
We simulate 400 HH neurons organized on a 20 × 20 grid with their probability of connectivity decaying exponentially with distance, capturing the exponential decrease in synaptic connectivity with cell distance as was done by McLaughlin et al. (2000). Then, the probability of connecting to a neuron on the grid a distance l in the horizontal direction and m in the vertical direction is given by where r = 20. The external drive to both networks is drawn from a Poisson process with rate ν and strength f Q , for the upstream, IAF network, and ν back and f back for the background drive to the downstream, HH network. This models input from neurons external to the patch of neurons that is explicitly being modeled. In addition, a fraction of the downstream population (20% of inhibitory cells and 30% of excitatory cells) receives an additional external drive modeling input from neurons external to the network that are firing in response to a sensory stimulus, mimicking the increased activity of a subset of neurons during a sensory experience. This additional sensory drive is modeled using a Poisson process with rate ν sens and strength f sens . Note that all figures were generated using 30 realizations of simulation runs for 5 s of simulated time in each trial. The inhibitory neurons are coupled with a gap junction with a coupling probability of 60% and neighboring excitatory pyramidal cells are coupled with a probability of 5%, reflecting experimental findings. One pair of electrotonically-coupled neurons, denoted the network-driven electrotonic pair (NDEP), is selected from the set of all electrotonically-coupled pairs in the network to be driven by a subset of neurons from the IAF network; the NDEP is denoted by the green and black triangles in Figure 1A. The IAF network is simulated prior to the HH model and input spikes are received by each neuron in the NDEP with strength f sens . The input neurons are chosen randomly from the IAF network such that ∼9-10 IAF neurons project to one neuron in the NDEP (the average firing rate received by each neuron in the NDEP is ∼100 Hz), with no IAF neuron projecting to both cells in the NDEP. The spike times of these input IAF neurons, shown in Figure 1B together with an example of the resulting dynamics of the downstream network in Figure 1C, are modeled as received by the NDEP in the same way as external spikes to network neurons (i.e., through a delta function in the G 4 variable at the input neuron spike times).
We postulate that synapses between the NDEP and the inhibitory cells are stronger than the average synapse between other cells in the network. The NDEP excites the gapjunction connected inhibitory cells to fire synchronous action potentials, which allows the synapses between the NDEP and the interneuron population to be preferentially enhanced through synaptic plasticity and classic learning rules (Feldman, 2012). As a result, we include a 10-fold increase in the synaptic strength from the NDEP to the interneurons and a reciprocal 3-fold increase in synaptic strength from the interneurons back to the NDEP, see Table 2 for normal synaptic strength values. This structure allows for a clear understanding of the effect of EC between a pair of pyramidal cells on network dynamics as discussed below.
The following sections describe the methods used to quantify the synchrony of the upstream IAF neurons' spike times, as well as the synchrony exhibited by the downstream network.

Input Synchrony
Input synchrony is a measure of the amount of synchronous input received by each neuron in the NDEP from the IAF upstream network. We determine it as follows. First, the instantaneous firing rate in each trial is calculated by binning time into 2 ms-sized bins and counting the number of spikes from the IAF network that occur in each bin (Shinomoto, 2010). This count is converted into an average firing rate by dividing by the Frontiers in Computational Neuroscience | www.frontiersin.org total number of neurons and the length of the time bin. Then, a smooth curve is created using a moving-window average with a smoothing width of 5 ms. We determine the value of the input synchrony by counting the total number of times in all trials that the smoothed instantaneous firing rate crosses a threshold of 35 Hz, i.e., the number of times the spikes from many neurons in the network fall within a short window of time. The threshold is chosen as the maximum of the smoothed instantaneous firing rate for a completely asynchronous network, yielding an inputsynchrony value close to zero for the asynchronous case, and guaranteeing that the measure will increase with an increase in synchrony. Note that this threshold is kept constant throughout all simulations and the results are not sensitive to moderate changes in this threshold. In this study, we adjust the input synchrony by varying the statistical properties of the Poisson external drive to the IAF network from fluctuation-dominated to mean-dominated, see Figure 2 for example networks with different input synchrony values. The input synchrony values are computed per trial. Clearly, the input synchrony measure correlates well with the coefficient of variation (CV) of the interspike interval, a measure for the regularity of the spiking statistics, see Figure 2C. However, we note that CV values give no indication of synchrony, or alignment of spikes in time, so we use the synchrony measure throughout the remainder of this work.

van Rossum Distance
The van Rossum method determines a distance, in time, between two pairs of spike trains by adding exponential tails to each spike time for each neuron and integrating the square of the difference between these exponential-tailed spike times (van Rossum, 2001). It is calculated as follows: First, the spike times of two neurons, x and y, are denoted as {t x i } and {t y j }, respectively. Then, exponential tails with time constant τ c are added to each spike time, t x i , to create a smeared-spike function, x(t), with the following equation: where H(t) is the Heaviside function. The smeared-spike function for the second neuron, y(t), is defined analogously. Finally, the van Rossum distance between the spike trains is defined as the integral of the square of the difference, The result is that D 2 is large for spike trains that do not align in time (are unsynchronized), and small for those that do (are synchronized). In this work, the time constant for the exponential tail, τ c , is set to 5 ms as was done in previous modeling studies ( Barranca et al., 2014) and is recommended for detecting synchrony on short time scales. Note that large values of τ C would yield an approximate difference in the total number of spikes elicited by the two neurons instead of the spike times.

Spike-Difference (SD) Method
We developed the spike-difference (SD) method to detect small changes in synchrony. Due to the gap-junction connectivity among the interneurons, the network exhibits synchrony for nearly all types of external drive, as well as in the presence or absence of a rare number of electrotonic pairs. The intuition is as follows: Once the network is sufficiently excited to induce firing, the interneurons fire synchronously due to their extensive gap-junction connections, causing a network synchronous event (NSE). Specifically, at the time of the NSE, the excitatory neurons are firing sufficiently to excite the interneuron population to threshold, creating a wave of inhibition that suppresses the activity of the entire population for a brief period of time succeeding the NSE. This can be seen in Figure 1B (right) by noting the lack of network spikes immediately after each NSE, where the NSE is indicated by the wave of inhibitory neurons spiking within a short window of time.
The SD method measures, for each NSE, the amount of synchrony exhibited by all neurons that participate in the event. The procedure is as follows. The time of the NSE is found by recording the times at which the magnitude of the average voltage of the interneuron population increases by 70 mV. This magnitude is kept constant throughout the simulations and the results do not change with moderate changes in this choice of magnitude, see Figure 3A for a raster plot of an example realization of the downstream network and the corresponding  Table 1 for values.
average voltage, with the threshold plotted as a purple line. Then, within a window of time around the NSE time, e.g., ±20 ms, the time-difference between each spike and the NSE time is measured for each neuron in each trial that spiked within this time window. The result is a distribution of time-differences in all trials, centered near 0 ms, the width, i.e., standard deviation, of which determines the average amount of synchrony the network exhibited during each NSE, see Figure 3B.

EC Enhances Synchrony Between a Pyramidal-Cell Pair
We begin by considering the effect of EC between one pair of neurons that is embedded in the downstream network and receives input from a subset of upstream IAF neurons exhibiting varying amounts of synchrony. First, we note that the two pyramidal cells in the NDEP have synchronized activity due to their EC, which strongly couples the voltages of the two neurons, resulting in almost identical spike trains. As an example, Figure 4A shows the spike times for the NDEP with experimentally-measured electrotonic conductance, compared with the spike times of the same two neurons with EC turned off, i.e., the junctional conductance is set to zero. Notice that the spike trains of the NDEP with EC turned on are almost identical, firing almost every spike together, whereas the synchronized firing between neurons in the NDEP is largely reduced with EC turned off.
We quantify the amount of synchrony exhibited by the spike trains in both cases (EC on and off) using the van Rossum distance, see Figure 4B. As expected from the near-identical spike trains, the distance between the spike times for the NDEP with EC turned on is significantly smaller than for the NDEP with EC turned off, for all input-synchrony values. Importantly, note that the structure and parameter values of the network with EC turned off are identical to the case with EC turned on, with the only difference between the two cases residing in the conductance of the EC between the two neurons in the one NDEP. The implication of this significant decrease is that downstream neurons will receive two synchronized excitatory spikes instead of two asynchronous ones, suggesting a role for the NDEP in organizing network activity around its own synchronized firing pattern.

One NDEP Increases Precision of Network Spike Timing
Due to the gap junctions that ubiquitously couple the network of inhibitory neurons, input from the excitatory population (whether the input is from neurons connected by a chemical synapse or by EC) creates synchronized firing events called NSEs (as described previously in the Methods section). The synchronization of neuronal spiking is often observed in response to sensory stimuli, with the timing of these synchronized events thought to be important for downstream communication (Azouz and Gray, 2003). The number of NSEs, which one can interpret as a network firing rate, is a relevant quantity of interest because it can also provide a measure for the amount of information transmitted from the model patch to a further downstream network. Figure 5A shows that the network containing an NDEP with EC turned on generates a greater number of NSEs than a network containing that same NDEP with EC turned off. Note that the NDEP with EC turned off still has the increased synaptic connection to and from the inhibitory network population, but that these two neurons exhibit significantly less synchronized spiking than two neurons with EC turned on, and so the inhibitory population does not often receive enough synchronized spikes to elicit an NSE in response to the firing of the NDEP. This means that the NDEP with EC turned off requires more input synchrony to elicit the necessary synchronized spiking to excite the inhibitory population and create an NSE. The NDEP with EC turned on, however, requires a less synchronous input to generate the same number of NSEs than an NDEP with EC turned off and thus we observe that the non-zero EC creates an NDEP that reduces temporal jittering of spikes and enhances the associations between spike events.
Firing rate is not the only method by which neuronal networks communicate and encode information. Rather, some neurons respond maximally to changes in spike timing on short timescales, i.e., tight synchrony (Singer, 1993). Figure 5B shows the SD measure using all neurons in the network (a measure of the tightness of the synchrony of each NSE, see the Methods section for the algorithm to compute the SD measure) to show that the NDEP with non-zero conductance creates NSEs that are tightly synchronized. The network in which the NDEP has zero conductance (EC is turned off) still generates NSEs; however, these events are less tightly synchronized than when EC is turned on, resulting in less communication to downstream areas. We also reinforce the result from the van Rossum distance (that the two neurons in the NDEP are tightly synchronized, recall Figure 4B), and validate the SD measure, by calculating the SD measure using only the spike times of the NDEP, see Figure 5C. Note that trend follows similarly as in the van Rossum distance, and we observe that the tight synchrony induced by the NDEP may aid in coincidence detection achieved by neurons in a downstream network.
These results together demonstrate that the NDEP with EC turned on elicits more NSEs than this same NDEP with EC turned off and, further, that these NSEs are composed of spike times that are closer together in time, i.e., more synchronized. However, notice that, in all three measures, EC turned on or off does not influence the behavior of the network as the input synchrony changes (see Figure 5 and notice that the lines appear nearly parallel). Therefore, in Figure 6, we quantify by averaging over all input-synchrony values the change in the van Rossum distance between the spike times of the neurons in the NDEP (Figure 4B), the synchrony (SD measure) of the neurons in the  NDEP for each NSE (Figure 5C), the total number of NSEs (Figure 5A), and the network synchrony (SD measure) of each NSE (Figure 5B) for the case when EC is turned on and off. These results combine to support the hypothesis that one NDEP with non-zero conductance can act to reduce the noise in an incoming stimulus (i.e., create a more synchronous event in response to less synchronous input) and increase the precision of network spike timing (i.e., tighten these synchronous events).

Two NDEPs Enhance Discriminability of Incoming Stimuli
Since Wang et al. (2010) identified several pairs of electrotonically-coupled pyramidal cells in the cortex, we next analyze a network in which an additional NDEP receives input from an upstream IAF network. Note that this additional NDEP is chosen randomly from the set of all electrotonic pairs in the network and has the same increase in synaptic strength to and from the interneuron network as in the one-NDEP case (see Methods section and Figure 1A). In addition, we note that the input to each NDEP can originate from disjoint sets of neurons in the same IAF network (SIAF) or from sets of neurons in two different IAF networks (DIAF) that have identical parameters (statistics), but differ in their incoming spike times (realizations of Poisson spike train). We investigate the dynamics of the downstream network in response to both types of IAF upstream input.
For the case in which two NDEPs receive DIAF input, we observe that the two NDEPs compete with one another through the inhibition evoked by each NSE, resulting in network behavior that is more variable than the case in which just one NDEP is present. Figure 7 demonstrates this competition by comparing the network containing two NDEPs (Figure 7A, top) to the network in which just one pair is chosen as the NDEP (Figure 7A, bottom), with the chosen pair indicated by the symbol "star" or "diamond." In comparing these networks, we determine how the behavior of each NDEP changes when an additional NDEP is present in the network. Notice that, in the network containing two NDEPs, there are several cases in which one NDEP fires, induces an NSE, and the wave of activity from the NSE inhibits the second NDEP from firing. This effect is shown in Figure 7B by the green circles in the bottom two panels, indicating those spikes that have been deleted by the competition between the NDEPs in the network receiving DIAF input.
The dynamics of a network containing two NDEPs differs from that containing one NDEP in several ways. First, we investigate how the average number of NSEs varies in each regime (i.e., EC turned off, one NDEP, two NDEPs receiving SIAF input, and two NDEPs receiving DIAF input). Figure 8A shows that a network containing two NDEPs generates a larger number of NSEs than a network with no electrotonic coupling (EC off) and a network with only one NDEP (EC on). Additionally, note that the number of NSEs generated by the two-NDEP network receiving DIAF input (purple bar) is smaller than the number generated by a hypothetical network in which there are no interactions between NDEPs (orange bar), calculated by doubling the number of NSEs in the one-NDEP case.   Figure 8A also demonstrates that the average number of NSEs generated by the two-NDEP network receiving DIAF input is larger than the number that is generated by the network in which the two NDEPs receive SIAF input (green bar). Figure 8B shows that this difference in the number of NSEs generated by each two-NDEP network becomes prominent when the input synchrony values become large (>25). The divergence of these two curves can be explained as follows. When two NDEPs are driven by the same IAF network, the spike times received by each NDEP become closer together as the input network synchronizes. As a result, the two NDEPs receive spikes at sufficiently close times so that their spike timing will be sufficiently close to elicit one NSE in response. In contrast, when two NDEPs are driven by different IAF networks, the spike times of the input neurons are sufficiently far apart so that their spike timing will also be sufficiently far apart to potentially elicit two distinct NSEs. This not only allows for the generation of more NSEs, but also creates the opportunity for the inhibition from one NSE to alter the spike timing of the other NDEP, as demonstrated in Figure 7. determining an excitatory event drawn as a black line and times at which excitatory events occurs indicated by gray-filled squares. The input IAF networks were driven with the following parameters for each regime: ν = 5, 000 Hz for EC off and one NDEP, ν = 4, 000 Hz for both two-NDEP networks.
Networks of neurons may communicate through changes in firing rate and synchrony; however, an important component of both types of communication is the pattern of this synchronized firing activity. The neurons that convey this information are the excitatory, or projection, neurons, whose axons innervate many networks. The activity of the excitatory neurons is modulated by the inhibition generated by each NSE, allowing for a variety of network firing patterns in response to the timing of the NSEs. Due to the competition between the NDEPs in the network containing two NDEPs receiving DIAF input, we observe a firing pattern that is much more variable than a network containing just one NDEP (see raster plots in Figure 7A for example). With more variability in the firing pattern of the network output, more information (more patterns of activity) can be transmitted to downstream areas.
To begin to quantify this variability, we examine the behavior of the excitatory-neuron population by looking at times at which a majority of the excitatory population is active, or the times at which the smoothed instantaneous firing rate (see Input Synchrony in Methods section for the algorithm to compute a smoothed instantaneous firing rate) of the excitatory population crosses a threshold of 1 Hz. For an example of the behavior of the excitatory neurons, for each regime, for one realization, see Figure 9.
Specifically, we are interested in the pattern of the timing of each of these excitatory events, see gray filled squares in Figure 9.
To observe the variability in the firing pattern, we compute the time difference between successive excitatory events, or interevent interval, for each regime. These time differences are then binned in a histogram with a bin size of 20 ms, see Figures 10A,B for all regimes.
FIGURE 10 | A network containing two NDEPs has greater variability in its firing pattern. (A) A histogram of the time difference between excitatory events, or inter-event interval, with EC turned off (pink), one NDEP (blue), and two NDEPs receiving DIAF input (purple). (B) A histogram of the time difference between excitatory events comparing the networks in which two NDEPs receive SIAF input (green) and DIAF input (purple). The coefficient of variation (CV), shown in the legend, is computed as the ratio of the standard deviation to the mean of the inter-event interval data. (C) The mean firing rate of each network averaged over 30 realizations (squares) together with the standard deviation (line). These results are for input synchrony values of 35 for both two-NDEP networks and 45 for the one-NDEP and EC-off networks, chosen to maintain a relatively constant firing rate across regimes, as shown in (C).
Notice that the distribution of the inter-event interval for the network containing two NDEPs (purple) is broader than for one (blue) and for the case when EC is turned off (pink), as shown in Figure 10A. The amount of variability in the activity of the network can be quantified using the coefficient of variation (CV), or the ratio of the standard deviation to the mean, with higher values indicating more irregularity in the time differences, see legend of Figure 10A. As expected, the CV for the case with two NDEPs receiving DIAF input is significantly larger than the CV for the other regimes. Note that the inter-event interval data for the two-NDEP case with DIAF input has a higher CV value (more variability) than the network receiving SIAF input due to the competition between the NDEPs (see Figure 10B), as discussed previously.
To confirm that the difference in the CV of the inter-event interval is not simply an effect of difference in firing rates, we show that the average firing rate of the neurons in each regime are similar, see Figure 10C. It is clear from the broad inter-event interval histogram, as well as the high CV value, that a network containing two NDEPs has the ability to output a wide variety of firing patterns, leading to a network that has the capacity to code for many different stimuli.

DISCUSSION
We have developed a model patch of a cortical network, organized in a way that elucidates a possible biological function for a rare number of electrotonic pairs, whose properties agree with experimental results. Measurements of EC between pairs of pyramidal cells in the adult cortex have been elusive to many experimentalists, having only recently been detected and their properties measured (Wang et al., 2010). Further investigations of our work (Crodelle et al., Manuscript in preparation) show that under normal network conditions, i.e., no strengthened synapses, a network with a rare number of electrotonic pairs is practically indistinguishable from one without; the observed effects are local and small, on the order of magnitude of standard network fluctuations. The rarity of these proposed ECs suggests that their role might be in a more specialized network, one in which the connection to the gap-junctionconnected interneuron population is enhanced. We have used such a network structure to show that EC between pairs of pyramidal cells could enhance the encoding of external stimuli. Specifically, we are interested in understanding how the NDEP (with EC turned on and off) alters the cortical network response to input that varies from unsynchronized to synchronized, with additional attention paid to analyzing how the varied HH-network output may be interpreted by further downstream areas.
In a network where synchrony is dominated by the gapjunction-connected interneuron population, we have shown that the addition of one NDEP generates many network-wide synchronous events by exciting the interneuron population with tightly synchronized excitatory spikes. These NSEs produce waves of strong inhibition, suppressing the activity of the network and keeping the firing rate of the network relatively low (sparse coding), while the tight synchrony of the NDEP keeps the timing of the NSEs very precise (spike-timing coding). This contrasts the network in which EC is turned off and an NSE can be induced only when a large percentage of excitatory neurons in the network are activated, resulting in high network firing rates and imprecise spike timing.
Coincidence detection in the cortex is an important property of pyramidal cells, in that it enables an efficient response to synchronized sensory input and organizes network activity (Azouz and Gray, 2000;Spruston, 2008;Shai et al., 2015). The tightly synchronized events elicited by the NDEP with EC turned on results in a network that sends very precisely-timed spikes to further downstream networks, aiding in coincidence detection capabilities of downstream areas. Our model predicts that networks in which EC between pairs of pyramidal cells is blocked will exhibit fewer network-wide synchronous events and have less precise spike-time correlation between incoming spikes from upstream networks, those elicited in downstream networks. From our preliminary numerical results, we also observe that the NDEP cannot be effectively replaced by just one excitatory neuron. This is because the spiking statistics of a pair of neurons with EC differs from those without EC in that pair of neurons, i.e., the neuron pair with EC elicits more regularly-timed spikes (lower CV values) and higher firing rates (smaller interspike intervals), when compared to that without EC.
In addition, studies have shown that the activity of interconnected inhibitory interneurons can modulate and influence the activity of the projection cells (Traub et al., 2001;Hjorth et al., 2009). Our investigation of a network with an additional NDEP shows that the two electrotonic pairs can elicit NSEs at their respective spike times, with a non-trivial interaction between the spike timing of each NDEP through inhibition from the interneuron population. The resulting network shows greater coding capacity than a network containing just one NDEP, as well as a network with EC turned off. Though we have only tested two NDEPs in this small-network case, our results predict that the coding capability for many NDEPs might be yet greater. The NSEs elicited by many NDEPs may create even more interactions between the spike timing of the NDEPs, yielding a larger variety of network firing patterns and further increasing the coding ability of the network. Although the size of our model network considered here is relatively small, the obtained results about the role for pair-wise EC in transmitting information across cortex through interactions with the vast gap-junction coupled interneuron populations are expected to extend to large-size neuronal networks.

DATA AVAILABILITY
All datasets analyzed for this study are included in the manuscript and/or the supplementary files. Computer codes and raw data are available upon request.

AUTHOR CONTRIBUTIONS
DC and GK contributed conception and design of the project. JC wrote the code, ran the simulations, and analyzed the results. DZ contributed to the interpretation of the results and generating new ideas. JC wrote the first draft of the manuscript. DZ, GK, and JC contributed to manuscript revision, read, and approved the submitted version.