Soma-axon coupling configurations that enhance neuronal coincidence detection

Coincidence detector neurons transmit timing information by responding preferentially to concurrent synaptic inputs. Principal cells of the medial superior olive (MSO) in the mammalian auditory brainstem are superb coincidence detectors. They encode sound source location with high temporal precision, distinguishing submillisecond timing differences among inputs. We investigate computationally how dynamic coupling between the input region (soma and dendrite) and the spike-generating output region (axon and axon initial segment) can enhance coincidence detection in MSO neurons. To do this, we formulate a two-compartment neuron model and characterize extensively coincidence detection sensitivity throughout a parameter space of coupling configurations. We focus on the interaction between coupling configuration and two currents that provide dynamic, voltage-gated, negative feedback in subthreshold voltage range: sodium current with rapid inactivation and low-threshold potassium current, IKLT. These currents reduce synaptic summation and can prevent spike generation unless inputs arrive with near simultaneity. We show that strong soma-to-axon coupling promotes the negative feedback effects of sodium inactivation and is, therefore, advantageous for coincidence detection. Furthermore, the feedforward combination of strong soma-to-axon coupling and weak axon-to-soma coupling enables spikes to be generated efficiently (few sodium channels needed) and with rapid recovery that enhances high-frequency coincidence detection. These observations detail the functional benefit of the strongly feedforward configuration that has been observed in physiological studies of MSO neurons. We find that IKLT further enhances coincidence detection sensitivity, but with effects that depend on coupling configuration. For instance, in models with weak soma-to-axon and weak axon-to-soma coupling, IKLT in the axon enhances coincidence detection more effectively than IKLT in the soma. By using a minimal model of soma-to-axon coupling, we connect structure, dynamics, and computation. Although we consider the particular case of MSO coincidence detectors, our method for creating and exploring a parameter space of two-compartment models can be applied to other neurons.


Introduction
Neurons that spike selectively to multiple subthreshold inputs that arrive within brief time windows are coincidence detectors. Coincidence detection is a fundamental neural computation that allows the brain to extract information from the temporal patterns of synaptic inputs. In the cortex, neurons have biophysical specializations compatible with coincidence detection [1][2][3][4][5][6][7], but some have questioned whether temporally-precise computations are possible in cortex due to highly variable neural activity therein [8,9]. In the early auditory pathway, the existence of coincidence detector neurons and their functional importance are widely valued [10][11][12]. Principal cells of the medial superior olive (MSO) in the mammalian auditory brainstem are a canonical example: they receive inputs originating from both ears [13,14] and are sensitive to microsecond-scale differences in the timing of arriving inputs [15][16][17][18]. These coincidence detector neurons are critical for sound-source localization [19] and likely play important roles in other aspects of binaural (two-eared) hearing such as sensitivity to interaural correlation [20,21] Temporally-precise neural coincidence detection requires specialized neural dynamics and circuitry. Coincidence detector neurons should have fast membrane dynamics with timescales of integration shorter than the intervals between volleys of synaptic inputs [2,22]. Inputs to coincidence detectors should also be brief and well-timed to precisely convey timing information. The requirements for effective coincidence detection in the auditory system are exceptionally stringent because auditory neurons must process inputs with temporal information at kilohertz-scale and higher [23][24][25]. Auditory brainstem circuitry is equipped with a suite of specializations to promote coincidence detection [26]. Afferent inputs to MSO cells are reliable and temporally-precise [27,28], dendritic processing in MSO further enhances coincidence detection [25,[29][30][31], and voltage-gated currents that are partially active near resting voltage make MSO cells extremely fast and precise processors [25,32]. Voltage-gated currents are also sources of dynamic negative feedback that contribute to the remarkable coincidence detection capabilities of these neurons. In MSO neurons, activation of low-threshold potassium (KLT) current and inactivation of sodium current are two identified sources of dynamic negative feedback [33][34][35][36]. In response to, say, a pair of brief excitatory inputs, these feedback mechanisms will transiently raise the spiking threshold after the first input, and thereby reduce the chance that the neuron will spike in response to the second input unless the inputs arrive nearly synchronously (within a coincidence detection time window).
We investigate the extent to which a structural specialization-namely, the coupling between the input region (soma and dendrite) and the output region (axon and axon initial segment)-can further optimize coincidence detection sensitivity. To do this, we develop a two-compartment neuron model as a minimal description of input-output coupling and systematically explore the effects of coupling configuration on coincidence detection sensitivity.
The two-compartment formulation is motivated by observations that spike generation in MSO likely occurs in the axon or axon initial segment [37,38]. Furthermore, sodium channels in the soma are inactivated near resting potentials [39] and spikes are small and graded in the soma [37], suggesting the soma does not participate in spike generation. Indeed, an absence (or small amount) of sodium current in the soma appears as a general design principle for temporally-precise auditory neurons [40]. Studies of coincidence detector cells in the avian auditory brainstem have shown that a passive soma can enhance coincidence detection [41] and that the distribution of voltage-gated channels in the axon initial segment undergoes activitydependent modulation [42,43] to improve coincidence detection, perhaps in an optimal manner [44]. A two-compartment formulation neglects the helpful contributions of dendritic processing to coincidence detection, but the role of dendrites has been considered in detail in previous studies [25,[29][30][31].
In this study, we systematically relate forward coupling (soma-to-axon) and backward coupling (axon-to-soma) strengths to model parameters. We explore this parameter space and find the coupling configurations that enhance coincidence detection sensitivity. Specifically, we identify strong soma-to-axon coupling as a natural configuration for neural coincidence detection because it engages sodium inactivation as a mechanism that transiently increases spike threshold on the time-scale of synaptic inputs and prevents firing to inputs that do not arrive concurrently. Moreover, the combination of strong soma-to-axon and weak axon-tosoma coupling generate spikes more efficiently (requires fewer sodium channels) and with shorter refractory periods than other models. This feedforward configuration enhances highfrequency coincidence detection and represents distinct advantages over one-compartment point neuron models that cannot exhibit this asymmetric coupling configuration. We observe that KLT current provides additional benefits for coincidence detection sensitivity, but these benefits depend on coupling configuration and where KLT current is located in the two-compartment structure. For instance, coincidence detection sensitivity in neurons with weak soma-to-axon coupling can be substantially improved if KLT current is co-localized with spike-generating currents.
We select passive properties to match known physiological characteristics of MSO neurons, so our observations apply directly to those canonical coincidence detectors. Nonetheless, our method for systematically exploring the parameter space of coupling configurations can be applied to study the relationships between structure, dynamics, and computation in other neurons that are well-described by a two-compartment idealization [24,[45][46][47][48][49][50]. In particular, a useful aspect of our work is that we show how to explicitly construct two-compartment models that satisfy the constraint of having nearly identical passive dynamics in the input compartment.
After a description of our model and method, we present our results as follows. First, we test coincidence detection sensitivity of the two-compartment model using simulated synaptic inputs. We then use simpler inputs (direct current) and phase plane analysis to explain the effects of coupling configuration on coincidence detection sensitivity. These main results concern interactions between cell structure and spiking dynamics. It is known, though, that KLT current is prominent in MSO neurons at subthreshold voltage levels and enhances their coincidence detection sensitivity. To complete our study, therefore, we repeat our study of coincidence detection sensitivity with simulated synaptic inputs, but this time with dynamic (voltage-gated) KLT current included in the model.

Two-compartment model parameterized by coupling strength
We construct and analyze a minimal description of a neuron that separates the input region (soma and dendrites) from the spike generating region (axon and axon initial segment) of a cell. This two-compartment model [45] has the form: The dynamical variables V i (i = 1, 2) describe the membrane potential in each compartment. Passive parameters in the model are membrane capacitance per area (C m ), axial conductance (g c ), reversal potential of leak current (E lk ), compartment surface area (A i ), and membrane leak conductance density (G i ). Parameters subscripted with i can take different values in the two compartments (i = 1, 2). To simplify notation, we will often omit the explicit reference to membrane area and instead use the notation c i = A i C m and g i = A i G i for i = 1, 2. The first compartment (abbreviation: Cpt1) receives input current (I in ) and the second compartment (abbreviation: Cpt2) is the site of spike-generating sodium current (I Na ). In some simulations we also include dynamic (voltage-gated) low-threshold potassium current (I KLT ). These currents are described in more detail below.
We use standard neurophysiological measures of passive activity in the soma to determine some parameters, and vary other parameters to create a family of two-compartment models distributed in a two-dimensional parameter space. We select model parameters so that, regardless of coupling configuration, the passive dynamics in Cpt1 are nearly identical regardless of the strength of coupling between compartments. This novel formulation allows us to meaningfully and systematically probe the dynamics of the model. The two parameters that define coupling configuration are introduced below. They describe strength of forward coupling (Cpt1 to Cpt2) and backward coupling (Cpt2 to Cpt1).
The properties we match to experimental measurements include resting potential in the soma (V rest ), input resistance for input to the soma (R in ), and exponential time constant (τ exp ) with which soma voltage returns to rest following a brief perturbation. We use the following values based on in vitro measurements of gerbil MSO neurons [32]: V rest = −58 mV, R in = 8.5 MO, and τ exp = 340 μs. We first match these properties using a model with passive dynamics (by setting I Na and I KLT to zero). After identifying parameter relations that satisfy these constraints, we discuss how to introduce sodium and KLT currents.
In the passive model, the resting potential is identical to the reversal potential and we have V rest = E lk . We now determine the remaining parameters based on the values of R in and τ exp . It is convenient to rescale the voltage equations by g i + g c and to introduce terms that represent the deviation of voltage from rest: U i = V i − E lk (for i = 1, 2). This yields the following equations for passive and subthreshold dynamics (I Na and I KLT removed, for now): The rescaled input current is denoted J in = I in /(g 1 + g c ). The time constants τ i = c i /(g i + g c ) describe the passive dynamics of the i th compartment (i = 1, 2) when the other compartment is held at its resting voltage.
We have also introduced in Eq 2 the two parameters that describe coupling strength. The forward coupling parameter is κ 1!2 . Formally, it is the ratio of voltages U 2 /U 1 at steady state in response to a constant current applied to Cpt1. Similarly, the backward coupling parameter is κ 2!1 . It is the ratio U 1 /U 2 at steady state in response to constant current applied to Cpt2. These quantities are attenuation factors that take values between zero (complete attenuation) and one (no attenuation). We find it more intuitive to refer to these constants as measures of coupling strength-values near zero represent weak coupling, values near one indicate strong coupling. We will refer throughout to κ 1!2 as the forward coupling parameter (soma-to-axon coupling) and κ 2!1 as the backward coupling parameter (axon-to-soma coupling). The relationship between coupling parameters and conductance parameters in Eq 1 are: Interestingly, the effect of one compartment on the other depends only on axial conductance and conductance in the target compartment. That is, the forward coupling parameter κ 1!2 depends on g c and g 2 but not on g 1 . Similarly, the backward coupling parameter κ 2!1 does not depend on g 2 . Moreover, axial conductance is the same for current flowing in either direction, so differences in forward and backward coupling are solely due to differences between leak conductance in the two compartments.
Next, we will show how to invert these equations to simply and uniquely define all passive model parameters for any combination of κ 1!2 and κ 2!1 . We only require prior knowledge of R in and τ exp (experimentally measurable parameters) and the assumption that the area of Cpt1 is much larger than the area of Cpt2. We denote the ratio of compartment areas as α = A 2 /A 1 . We will always use α = 0.01 in this study. This assumption is plausible for cells with input regions that are much larger than spike-generating regions, and is consistent with previous models of auditory coincidence detector neurons [24,38].
We find the axial conductance (g c ) by expressing it in terms of input resistance and the coupling coefficients. By setting U 0 1 ¼ U 0 2 ¼ 0 in Eq 2, we find the steady state relations U ss 1 ¼ k 2!1 U ss 2 À J in and U ss 2 ¼ k 1!2 U ss 1 . From these relations, and the steady-state input resistance (applying Ohm's Law), it follows that R in ¼ À U ss Solving for g c and after some substitutions we find The remaining parameters determine the passive dynamics of V 1 , so they depend on τ exp . In a one-compartment passive model, τ exp is identical to the membrane time constant and its value is c 1 /g 1 . In a two-compartment model, coupling between the compartments introduces a second time-scale that can influence the rate at which V 1 returns to rest after a brief perturbation (Rall's equalization time constant [51], and see also [52]). Some additional analysis is required, therefore, to relate τ exp to model parameters.
We invoke the assumption that the input region is much larger than the output region (A 1 � A 2 , or, equivalently, take 0 < α � 1) and observe that this can create a separation of time-scales in the passive dynamics of U 1 and U 2 . The ratio of time constants in the two compartments is τ 2 /τ 1 = (c 2 /c 1 )[(g 1 + g c )/(g 2 + g c )]. After some substitutions, and using the assumption that C m is identical in both compartments, we find that We restrict ourselves to coupling configurations for which κ 1!2 /κ 2!1 does not exceed ten, so that τ 1 is an order of magnitude larger than τ 2 (recall we use α = 0.01). In this scenario, we can segregate the passive dynamics into a slow variable (U 1 ) and a fast variable (U 2 ). The ratio of time-constants is a small parameter which we denote by � ¼ a k 1!2 k 2!1 . For � close to zero, we can make the approximation that U 2 evolves instantaneously (on the fast time-scale) to its U 1dependent steady-state value of U ss On the slow-time scale, U 2 takes this instantaneous value and the dynamics of U 1 are (to leading order in the small parameter �): In other words, in cases when τ 1 � τ 2 , the passive dynamics in Cpt1 are approximately linear and one can match the time-scale of U 1 to the experimentally-observed membrane decay time by setting τ 1 = τ exp (1 − κ 1!2 κ 2!1 ). This slaving of U 2 to U 1 is valid for describing passive subthreshold dynamics. If sodium current is included, then the dynamics are non-linear and spike-generation is possible. On the slow time-scale, spike-generation is represented by a discontinuous jump to a fixed point at higher values of U 1 and U 2 .
To summarize our method: we use values for three standard neurophysiological measures of (passive) soma dynamics (R in , τ exp and E lk ), we choose α (the ratio of surface areas A 2 /A 1 ) to be small (α = 0.01 in all simulations), and we let the two coupling constants define a twodimensional parameter space of soma-axon coupling. For any coupling configuration, we can then uniquely determine the passive parameters in Eq 1. The parameter relationships, described above, are: Using these parameter relations guarantees that passive dynamics in Cpt1 remain nearly identical as we explore neural dynamics and coincidence detection sensitivity in the twodimensional parameter space of coupling strengths.
We only consider coupling configurations in which forward coupling is stronger than backward coupling (κ 1!2 � κ 2!1 ). This corresponds to an assumption that voltage signals propagate forward from the soma to the axon with less attenuation than signals that backpropagate from the axon to the soma. This condition is appropriate for MSO neurons, since in vitro recordings show weak backpropagation of action potentials to the soma and dendrites [31,37]. Notice, from Eq 3, that the condition κ 1!2 � κ 2!1 is met if g 1 � g 2 . This is reasonable in our model because we assume A 1 � A 2 and do not expect that the density of leak conductance in each compartment (G 1 and G 2 ) would differ by orders of magnitude.

Low-threshold potassium model
A voltage-gated low-threshold potassium current (I KLT ) is prominent in MSO neurons and thought to improve coincidence detection sensitivity [53]. We model I KLT in the i th compartment with the equation where the reversal potential is E K = −106 mV. We include the second term so that the addition of KLT current does not alter the resting potential (I KLT,i = 0 when V i = V rest , for i = 1 or 2.). Equivalently, one could adjust E lk to counterbalance the amount of KLT current active at rest, and omit this correction term. The dynamics of the activation variable w i are as in [32] at a temperature of 35˚C: The inactivation variable is slow (time-scale of several hundred milliseconds), so we make the simplification that its value is fixed at the steady state z 1 (V rest ) where the steady state function is [32] Next we discuss how we include dynamic I KLT in the two-compartment model. We omit subscripts, for ease of presentation, but the same method applies for dynamic KLT current in either compartment. We first find the passive leak conductance in the relevant compartment using Eq 7. Call this g lk . We then reduce this conductance by some amount, typically 10%. In other words, we set the leak conductance in the relevant compartment (g i ) to 0.9g lk . Lastly, we set g KLT to the value that preserves the total conductance in the compartment at the resting potential. In some simulations, we leave the KLT conductance fixed at its resting value. We refer to this case as frozen KLT-the KLT current acts as a leak current and the subthreshold dynamics are the same as the original passive model. In other simulations, we allow KLT conductance to depend on voltage. We refer to this case as dynamic KLT. To include dynamic KLT in Cpt1, for example, we would choose g KLT so that it satisfies the equation g 1 þ g KLT w 4 1 ðV rest Þz 1 ðV rest Þ ¼ g lk and allow the KLT activation variable w to evolve according to Eq 9.

Sodium current model
The second compartment represents regions of the cell in which spikes are generated, presumably the axon initial segment or other excitable regions in the axon [38]. We use a reduced model of sodium current, adapted from earlier models of auditory brainstem neurons [32,54], to produce spikes: I Na ðV 2 ; hÞ ¼ g Na m 3 1 ðV 2 ÞhðV 2 À E Na Þ À g Na m 3 1 ðV rest Þh 1 ðV rest ÞðV rest À E Na Þ; where the sodium reversal potential is E Na = 55 mV. The second term is included so that Coupling configurations that enhance coincidence detection I Na = 0 when V 2 is at its resting value. Equivalently, one could adjust E lk to counterbalance the amount of sodium current active at rest, and omit this correction term. Setting I Na = 0 at rest simplifies analysis of the model and is appropriate for MSO neurons. Sodium current in MSO neurons is small at rest, with most sodium channels inactivated at the resting membrane potential [39]. We assume that activation of sodium is sufficiently fast to justify the approximation that the gating variable m instantaneously reaches its voltage-dependent equilibrium value m 1 (V 2 ) [34]. The gating variable h governs inactivation of the sodium current and has dynamics as in [32]: These are the same as in [54], but with temperature adjusted to 35˚C, and also note the resting membrane potential in our model is -58 mV as opposed to -65 mV in [54] Our primary objective is to determine the effects of coupling configuration on coincidence detection. Maximal conductance g Na determines excitability and spike threshold in the model neuron and thus also influences coincidence detection sensitivity [22]. Rather than setting g Na to an arbitrarily chosen value, we explore a range of g Na to determine the best possible coincidence detection sensitivity over this range of g Na , for each coupling configuration. We explain our method for choosing g Na values in more detail below.

Synaptic input model
We generate synaptic inputs to the two-compartment model using a model of the auditory periphery [55]. This model includes the effects of cochlear filtering and nonlinearities, inner hair cell activity, and synaptic transmission, and generates auditory nerve spike trains. As inputs to this model, we use sine waves that represent pure tone sounds. We perform simulations with frequencies ranging from 200 Hz to 700 Hz at a level of 70 dB. The neuron model receives two streams of auditory nerve inputs representing (conceptually) inputs from the two ears, see the schematic in Fig 1. The sine waves to the two ears are presented either with identical timing to generate coincident inputs, or with a time delay to simulate non-coincident inputs.
MSO coincidence detectors receive a small number of synaptic inputs [56], so we use the auditory nerve model to simulate five independent input sequences of spike times per ear. Each auditory nerve spike time creates an excitatory post-synaptic conductance (EPSG) described by a double exponential function [18]: Conductance g syn is in units of nanosiemens and t is milliseconds. These EPSGs are transformed into synaptic current (EPSCs) according to the equation where the reversal potential for the excitatory current is E syn = 0 mV. We set the constant scaling factor in the definition of g syn (t) (Eq 13) so that a single excitatory input depolarizes V 1 by roughly 6 mV, a value consistent with measurements of MSO neurons' responses to synaptic excitation in vitro [57]. Our goal is to identify and understand essential aspects of how structural and dynamical features of neurons affect coincidence detection. Although we do not incorporate a complete description of neural processing in the MSO-pathway, we view this setup to be adequate for probing coincidence detection in an MSO-like two-compartment model using quasi-realistic stimuli. Notably, we do not include spherical bushy cells in the cochlear nucleus that may enhance temporal precision of afferent inputs to MSO neurons [27], nor do we include inhibitory inputs that appear to modify time-difference tuning of MSO coincidence detectors [57][58][59][60][61][62][63], but see also [18].
In some simulations, as we will make clear in the context of the Results, we use more simplistic inputs such as steps or ramps of current injected directly into the input compartment to study response characteristics of the model. . Sine-wave inputs are delivered to a model of the auditory periphery [55]. We extract from this model simulated spiking responses of five auditory nerve fibers for each sine-wave. These simulated spike trains are used to create excitatory current (I in in this figure) that is delivered to Cpt1. B: Example synaptic conductance g syn for 500 Hz stimuli that are either in-phase (left panel) or out-of-phase (right panel). Gray line is the amplitude of a unitary conductance event. https://doi.org/10.1371/journal.pcbi.1006476.g001

Measure of coincidence detection sensitivity
We measure firing rate (spikes per second) generated by the two-compartment model in response to coincident inputs (identical sine wave stimuli to the two ears of the auditory nerve model) and non-coincident inputs. In some simulations we generate non-coincident inputs by using sine wave stimuli that were anti-phase to the two ears of the auditory nerve model. For example, for a 500 Hz stimuli, the two sine waves would have a 1 ms time difference. In this construction of non-coincident inputs, the time difference between the sine waves shortens with increasing frequency. To confirm that this dependence of time difference on frequency does not bias our results, we also perform simulations in which we generate non-coincident inputs by using sine wave stimuli to the auditory nerve model with a fixed time difference of 500 μs.
To measure coincidence detection sensitivity, we compute the difference in firing rates for responses to coincident and non-coincident inputs. We compute firing rates (spikes per second) by counting the number of spikes generated in trials of length 250 ms. We then calculate mean and standard error of firing rates from 100 repeated trials. Large-scale calculations to sweep over parameter space were performed using Matlab simulation code executed on computers managed by the Ohio Supercomputer Center. The ordinary differential equations defining the two-compartment model were solved numerically using the Matlab command ode15s (a variable-step, variable-method solver useful for stiff systems). Simulation code is available at https://github.com/jhgoldwyn/TwoCompartmentModel.
A good coincidence detector neuron would be one with a large difference in firing rates for these two conditions. Firing rate difference measures of coincidence detection sensitivity have been used in related studies [64,65]. Other measures have been considered, including Fisher information [38], width of time-difference tuning curves [62], and quality factor (similar to dprime) [22]. The right measurement of coincidence detection sensitivity remains, as these alternatives reveal, an open question (and one wrapped up in ongoing debates regarding the nature of the neural code for sound source location [66]).
One justification for comparing in-phase to out-of-phase firing rates is that it is relevant to a system that uses a two channel representation of auditory space in which sound location is represented by the difference in firing rates between two populations of cells tuned to distinct time-differences [67]. We suggest an additional perspective based on an analogy to signal classification theory and the receiver operating characteristic (ROC) [68].
Consider a coincidence detector neuron responding to a periodic (sine wave) stimulus. Each cycle of the stimulus evokes a volley of synaptic inputs that may or may not be temporally aligned with one another. The task of the coincidence detector neuron is to respond (generate a spike) if the synaptic inputs arrive within a brief time window and to not respond (not spike) if the synaptic inputs are dispersed in time. From this perspective, a coincidence detector neuron is an observer of its own synaptic inputs and it signals the presence of coincident inputs by generating a spike. Chance [69] has articulated a similar approach for measuring synaptic efficacy.
Extending the analogy, for each two-compartment model (parameterized by coupling configuration), we construct ROC curves by plotting hit rate (firing rate to coincident inputs) against false alarm rate (firing rate to non-coincident inputs) for varying values of the sodium conductance g Na . Sodium conductance controls the overall excitability of the model and operates as the threshold parameter in ROC analysis. To compare coincidence detection sensitivity across coupling configurations, we simulate the model for a range of g Na values and define coincidence detection sensitivity to be the maximum firing rate difference. In this way, we identify the g Na level for which the neuron, acting as an observer of its inputs, is the best possible coincidence detector. In other words, we identify the g Na value that maximizes hits (spikes generated in response to coincident inputs) while minimizing false alarms (spikes generated in response non-coincident inputs), for a given coupling configuration and stimulus. We illustration this calculation in Fig 2A (cartoon only, not simulation data). We provide S1 Fig. to illustrate this calculation with simulation data.
The g Na values that produced similar spiking activity for different coupling configurations could vary by orders of magnitude, so we required a reasonable way to determine the appropriate range of g Na values to use. For each coupling configuration, we set a reference value for g Na , denoted g ref Na , by finding the smallest g Na value for which a pair of coincident EPSGs could evoke a spike. The rationale for this definition of g ref Na is that, for any g Na value larger than g ref Na , the model neuron could possibly spike in response to two non-coincident EPSGs. For g Na significantly larger than g ref Na , the model neuron could even spike in response to a single EPSG. By restricting g Na to values near g ref Na  From the firing rate differences measured across this range of g Na values, we identify the maximum firing rate difference and use this as our measure of coincidence detection sensitivity. As a result of this process, we obtain different best g Na values depending on input frequency, coupling configuration, and KLT currents. An implication of this approach is that neurons should modulate sodium conductance based on stimulus parameters and physiological conditions. We do not pursue this idea here, but see [42,44] for relevant studies and S3 Fig. for supporting results showing changes of best g Na with stimulus frequency. Firing rate to coincident inputs and noncoincident inputs increase with g Na (cartoon, not actual data). We sweep across a range of g Na values and quantify coincidence detection sensitivity as the maximum firing rate difference across all g Na values used (inset). We draw an analogy to the signal receiver operating characteristic (ROC) curve: coincident firing rate is the hit rate, noncoincident firing rate is the false alarm rate, and g Na sets the detection threshold. This panel is for illustration only and does not portray actual simulation data. B: Reference values for sodium conductance (g ref Na ; the smallest value of g Na at which a pair of simultaneous EPSG events evokes a spike). The upper left region is empty (white color) in this and subsequent plots of the space of coupling parameters because we only consider models with forward coupling stronger than or equal to backward coupling (κ 1!2 � κ 2!1 , see text for discussion). For further analysis of the dynamics of the two-compartment model, we performed bifurcation analysis using the Auto feature in XPPAUT. XPPAUT is a software package for analyzing and simulating dynamical systems, and uses continuation methods to perform numerical bifurcation analysis [70].

Parameterization of a family of two-compartment models by coupling strength
As discussed in Materials and methods, the forward coupling parameter (soma-to-axon, κ 1!2 ) and backward coupling parameter (axon-to-soma, κ 2!1 ) characterize the passive dynamics of a two-compartment model (see Eq 3). By defining model parameters according to Eq 7, we construct a family of models with nearly identical passive dynamics in Cpt1. Recall we select model parameters based on reported properties of MSO neurons: steady state input resistance is 8.5 MO, the decay time constant is 340 μs, and the resting potential is −58 mV [32]. We also assume the capacitance per unit area is the same in both compartments and that the surface area of Cpt1 was 100 times larger than the surface area of Cpt2 (α = 0.01). Using a typical value of membrane capacitance C m = 0.9 pF/cm 2 [71], we find that capacitance in Cpt1 is C 1 = 40 pF. The area of Cpt1 is A 1 = 4444 μm 2 , which is a plausible value for the surface area of the soma and dendrite regions of an MSO neuron [31,72].
In order to maintain identical passive dynamics in Cpt1 across coupling configurations, the leak conductance (g 1 , g 2 ) and axial conductance (g c ) vary with the values of the coupling parameters as shown in Fig 3. Leak conductance in Cpt1 increases as forward coupling strength increases but decreases as backward coupling strength increases (Fig 3A). Leak conductance in Cpt2 exhibits the opposite dependence on coupling strength: g 2 decreases as forward coupling strength increases and increases as backward coupling strength increases ( Fig  3B). The axial conductance connecting the two compartments depends primarily on the strength of backward coupling-the contour lines in Fig   Parameter values for Cpt1 leak conductance g 1 (Panel A), Cpt2 leak conductance g 2 (Panel B), and axial conductance g c (Panel C). These parameter values are derived using Eq 7 and the assumptions that the decay time constant for V 1 is τ exp = 0.34 ms, input resistance in Cpt1 is R in = 8.5 MO, and the surface area of Cpt1 is one hundred times larger than the surface area of Cpt2. The upper left half of the parameter space is empty because we did not consider models for which forward coupling was weaker than backward coupling. Colored stars mark the locations of three configurations we examine in detail below: the weakly-coupled model (κ 1!2 = 0.3, κ 2!1 = 0.2; blue), the forward-coupled model (κ 1!2 = 0.8, κ 2!1 = 0.2; green), and the strongly-coupled model (κ 1!2 = 0.8, κ 2!1 = 0.7; red). https://doi.org/10.1371/journal.pcbi.1006476.g003 Coupling configurations that enhance coincidence detection panel of Fig 3 is empty because we only consider coupling configurations for which forward coupling is not weaker than backward coupling (κ 1!2 � κ 2!1 ).
To explore how coupling configuration modifies neural dynamics, we will often compare three models near the edges of the coupling parameter space. These are a weakly-coupled model (κ 1!2 = 0.3, κ 2!1 = 0.2), a strongly-coupled model (κ 1!2 = 0.8, κ 2!1 = 0.7), and a forward-coupled model (κ 1!2 = 0.8, κ 2!1 = 0.2). The locations of these three models in the coupling parameter space are shown as colored stars in Fig 3. We remark that complete coupling (κ 1!2 = κ 2!1 = 1) is equivalent to a one-compartment point neuron model because voltages in the two compartments are the same in this case.

Passive dynamics in first compartment are nearly identical across coupling configurations
Our parameterization method is designed to maintain the same voltage response in Cpt1 (V 1 ) regardless of the coupling configuration. In fact, due to the strong separation of time scales between the two compartments (recall Eq 5), the voltage in Cpt1 is governed approximately by linear dynamics with time constant τ exp (see Eq 6) and the voltage in Cpt2 is . These approximations are valid to leading order in the small parameter � ¼ a k 1!2 k 2!1 . We remind the reader that these calculations are performed in the case of passive dynamics-i.e. for a model without spike-generating sodium current (g Na = 0) and with frozen low-threshold potassium current (I KLT acts as a leak current, and in fact is equivalent to g KLT = 0, see Materials and methods).
Simulations of passive two-compartment models illustrate how the parameterization method results in models with nearly identical V 1 dynamics (Fig 4B). We use darker colors to show time-courses of voltage in Cpt1 (V 1 ) and lighter colors to show time-courses of voltage in Cpt2 (V 2 ). V 1 responses are nearly identical regardless of coupling configuration and attenuation of V 2 responses depends on κ 1!2 . Time-courses of V 1 and V 2 shown here are responses to 500 Hz coincident inputs. The three coupling configurations in this figure are indicated in the schematic in the top row (from left to right in Fig 4A: weakly-coupled, forward-coupled, and strongly-coupled, as defined previously).

Spiking dynamics depend on coupling configuration
We include spiking in the model by adding sodium current to Cpt2 (see Eq 11). As described in Materials and methods, we define g ref Na for each coupling configuration as the minimum level of g Na at which two coincident inputs evoke a spike. We find it helpful to normalize g Na to these reference values when comparing across coupling configurations. In Fig 4C, we show responses to 500 Hz coincident inputs with g Na set to g ref Na . This results in g Na = 6291 nS for the weakly-coupled model, g Na = 398 nS for the forward-coupled model, and g Na = 2003 nS for the strongly-coupled model. We do not include dynamic KLT current in these simulations.
These responses to identical inputs (same trains of synaptic inputs for each model, regardless of coupling configuration) reveal that spiking dynamics differ depending on coupling configuration. Coupling configuration dramatically changes the amplitude and shapes of spikesthe peaks of V 2 in the weakly-coupled model are near 40 mV whereas the peaks of V 2 in the forward-coupled and strongly-coupled models are near 0 mV (note the different axes). Spikes in the weakly-coupled model tend to be all-or-none, but the forward-coupled and stronglycoupled models can have more graded spike amplitudes. Not surprisingly, amplitudes of backpropagated spikes in Cpt1 depend on κ 2!1 . For models with weak backward coupling (the weakly-coupled and forward-coupled models), V 2 spikes cause small V 1 deflections that can appear similar to subthreshold V 1 activity. For the strongly-coupled model, in contrast, V 1 tracks V 2 more closely and shows larger amplitude backpropagated spikes in Cpt1.
Importantly, the number and timing of spikes also changes with coupling configuration. In these traces, the forward-coupled model has two more spikes than the weakly-coupled and strongly-coupled models (see the extra spikes at 10 ms and 20 ms for the forward-coupled model). This anticipates our main result: coupling configuration affects spike generation and can alter the sensitivity of neurons to coincident inputs. If we view this example simulation using the analogy to signal detection theory and ROC analysis, described previously, we can say that the forward-coupled model correctly identifies two more coincident events (has two more hits) than the weakly-coupled and strongly-coupled models.

Optimal coupling configuration for coincidence detection in twocompartment models with frozen KLT current
In a first set of simulations, we study the two-compartment model with passive subthreshold dynamics. The low-threshold potassium (KLT) current is frozen and included as part of the leak current (see Materials and methods). The only voltage-gated current in this set of simulations is the spike-generating sodium current in Cpt2. We quantify coincidence detection sensitivity by finding the maximum firing rate difference between coincident and noncoincident inputs for each coupling configuration (as described in Materials and methods). In Fig 5, we report results for three stimulus frequencies (from left to right: 300 Hz, 500 Hz, and 700 Hz). We construct non-coincident inputs in two ways: in Fig 5A we use out-of-phase sine wave inputs to the auditory nerve model. Time delays for out-of-phase inputs vary with frequency, so we also test coincidence detection sensitivity using sine wave inputs that are misaligned in time by a fixed 500 μs time difference in Fig 5B. Larger firing rate differences reflect better coincidence detection. In our analogy to a signal classification task, we say a large firing rate difference indicates a detector with a high hit rate in response to coincident inputs and a low false alarm rate in response to non-coincident events. We make two observations that we will explore in greater detail below. First, coincidence detection sensitivity improves with increases in forward coupling strength. That is to say, as one moves from left-to-right within each panel of Fig 5, the firing rate difference increases. Second, the combination of strong forward coupling (large κ 1!2 ) and weak backward coupling (small κ 2!1 ) enhances coincidence detection for high-frequency stimuli (notice Coincidence detection sensitivity measured as the maximum firing rate difference between responses to in-phase and out-of-phase stimuli. B: Coincidence detection sensitivity measured as the maximum firing rate difference between responses to inphase stimuli and those with a 500 μs time difference. In both rows: coincidence detection sensitivity is measured using three different stimulus frequencies (from left to right: 300 Hz, 500 Hz, 700 Hz). Color scheme in each panel is normalized so that the lowest color (dark blue) is the firing rate difference that is 65% of the highest value in that panel. In many panels these values are not reached and range of colors presented does not include these low blue colors. https://doi.org/10.1371/journal.pcbi.1006476.g005 Coupling configurations that enhance coincidence detection the large firing rate differences in the lower right corner of panels for 500 Hz and 700 Hz stimuli).
In Fig 6, we further detail how coincidence detection sensitivity changes with stimulus frequency. Strong forward coupling enhances coincidence detection for all frequencies above 200 Hz (the weakly-coupled model has the smallest firing rate difference across these frequencies). An advantage for the forward coupling configuration (combination of strong forward coupling and weak backward coupling) emerges for stimulus frequencies above about 400 Hz. These observations hold (qualitatively) regardless of whether we use out-of-phase stimuli for non-coincident inputs (Fig 6A) or time-delayed stimuli ( Fig 6B). We provide tuning curves in S4 Fig. for additional views of how coincidence detections sensitivity depends on stimulus frequency, g Na , and coupling configuration.
We emphasize that the comparisons we are making are across coupling configurations. We are not suggesting, in particular, that the non-monotonic trend in Fig 6A indicates an optimal frequency for coincidence detection. This non-monotonic trend is due in part to a ceiling effect at low frequencies (in-phase firing rates do not exceed 200 Hz for a 200 Hz stimulus), and also masks the fact that the proportion of coincident inputs that trigger a spike (the hit rate, in the signal detection analogy) decreases with increasing stimulus frequency.
Why does strong forward coupling improve coincidence detection? And why does the specific combination of strong forward and weak backward coupling (the forward-coupled model) enhance high-frequency coincidence detection? We provide explanations below. First, we will demonstrate that strong forward coupling endows the two-compartment model with two properties that are advantageous for neural coincidence detection: phasic firing and sensitivity to input slope. Second, we will show that the specific combination of strong forward coupling and weak backward coupling shortens the refractory period of the two-compartment model. Neurons with short refractory periods can faithfully and rapidly respond to highfrequency sequences of coincident inputs. Neurons with longer refractory periods are Coupling configurations that enhance coincidence detection disadvantaged when performing high-frequency coincidence detection. They may miss opportunities to generate spikes in response to coincident inputs, and thus their hit rate (when thinking of these neurons as signal detectors) may be depressed.
Strong forward coupling ensures phasic firing. Here we show that coupling configuration determines whether two-compartment models respond to steps of current injection with a single spike (phasic firing) or periodic spike trains (tonic firing). Phasic firing is advantageous for coincidence detection because it allows neurons to respond to rapid changes in stimuli-the concurrent arrival of several synaptic inputs, e.g.-and remain insensitive to slower changes in stimulus level that do not signal coincident inputs.
In these simulations we use direct current injection (not synaptic inputs) and frozen KLT current. See Fig 7A for example waveforms of injected current. In Fig 7B-7D, we show timecourses of voltage in Cpt2 to these step current inputs, with g Na set to g ref Na . As the amplitude of the applied current increases, the weakly-coupled model exhibits a transition from quiescence (no firing) to periodic firing (Fig 7B). The transition occurs via a Hopf bifurcation and is a The lowest values of g Na at which tonic firing is possible is g tonic Na . Models with strong forward coupling and weak backward coupling (lower right region in this panel) exhibit phasic firing over the largest range of (normalized) g Na values. tonic firing pattern. The forward-coupled and strongly-coupled models do not fire repetitively in these simulations. Instead, they transition from quiescence to a single spike elicited at the onset of the current injection (Fig 7C and 7D). This response pattern is defined as phasic firing. We remind readers familiar with the Rothman-Manis [54] model that, in addition to introducing a two-compartment structure, we have also made changes to parameter values to reflect MSO physiology (time constant, input resistance, and resting potential). Thus, our observation of phasic firing in the strongly-coupled model does not contradict the observation by [36] of tonic firing in a Rothman-Manis-like model with frozen KLTcurrent, even though a one-compartment model can be viewed as an extreme case of a two-compartment model with κ 1!2 = κ 2!1 = 1. The mechanism for phasic firing in our model with frozen KLT is the dynamic sodium inactivation gating variable h [33,36]. We discuss this in more detail below, when we analyze the dynamics of the system in the V 2 -h phase plane.
We explore a range of g Na values and trace out the boundary in parameter space between tonic and phasic firing. For each coupling configuration, we vary g Na and the amplitude of the steady current (I step ). We performed a two-parameter bifurcation analysis. The regions within the U-shaped curves in Fig 7E are parameter combinations (g Na and current levels) that produce repetitive firing.
The two-compartment model exhibits a tonic firing pattern if g Na is set to a sufficiently large value. We identify the smallest values of g Na at which repetitive firing to steady current can be observed (the lowest point on each U-shaped curve), and label these values g tonic Na . For example, the weakly-coupled model can exhibit tonic firing for values of g Na larger than 0:5g ref Na , roughly. The forward-coupled and strongly-coupled models are phasic unless g Na is larger than 3:3g ref Na or 2:7g ref Na , respectively. Values of the ratio g tonic Na =g ref Na for all coupling configurations are shown in Fig 7F. The value of this ratio increases with forward coupling strength. Phasic dynamics are closely associated with temporally-precise neural coincidence detection [34,35]. This result indicates that phasic dynamics in models with stronger forward coupling are robust, in the sense that this coupling configuration can allow neurons to maintain phasic dynamics over a larger range of g Na and input levels.
Strong forward coupling enhances slope sensitivity. Another property associated with coincidence detection is slope sensitivity. Slope-sensitive neurons are those for which the rate of increase (slope) of an input-not just the input amplitude-can determine whether the neuron generates a spike [34]. As we now demonstrate, models with strong forward coupling are responsive to inputs with fast-rising inputs, whereas models with weak forward coupling can respond to inputs with slow-rising slopes. This is a second indication that strong forward coupling benefits neural coincidence detectionin models with frozen KLT current.
We vary the slope of inputs using ramps of current, as shown in Fig 8A. Time-courses of V 2 In response to these three ramps are shown in Fig 8B-8D. We fix the maximum amplitude at 2500 pA for these inputs, and let the duration of the ramp vary with ramp slope.
These simulations demonstrate that the forward-coupled model is the most slope-sensitive, as it only generates an action potential in response to the ramp with the steepest slope. In contrast, the weakly-coupled model spikes in response to all three ramps, and the strongly-coupled model spikes in response to two out of the three ramps.
We measure slope sensitivity by finding the smallest g Na at which each model spikes in response to ramps of varying slopes (Fig 8E). We again fix the amplitude of the current ramp at 2500 pA. The stars in Fig 8E mark the parameter values used in Fig 8B-8D. For the weaklycoupled model, this measure of slope threshold is relatively constant as a function of ramp slope. In contrast, g Na at firing threshold decreases dramatically for the forward-coupled and strongly-coupled models. This confirms that the forward-coupled and strongly-coupled models are more sensitive to input slopes than the weakly-coupled model. By setting g Na appropriately, for example g Na � g ref Na , the forward-coupled and strongly-coupled models can be tuned to selectively spike in response to rapidly-rising inputs but not gradually-rising inputs.
To compare slope-sensitivity across all coupling configurations, we measured g Na at firing thresholds in response to slow and fast rising ramps (500 pA/ms and 1000 pA/ms) and defined Δg Na to be the ratio of these values. A large value of Δg Na indicates that the model neuron can be made to fire selectively to steeply-rising ramps only (a slope-sensitive neuron), whereas values of Δg Na near one indicate that the model neuron responds similarly to ramps regardless of input slope. By this measure, models with strong forward coupling are more slope-sensitive than models with weak forward coupling. Moreover, the forward-coupled configuration (lower right corner of Fig 8F) is optimal for slope-sensitivity. The strongly-coupled model responds with a spike in response to the two ramps with steeper rising slopes. E: Sodium threshold is the smallest value of g Na (plotted after normalization by g ref Na ) for which a ramp of a given slope elicits a spike. Threshold for the weakly-coupled model is relatively constant with ramp slope, but decreases with ramp slope for the forward-coupled and strongly-coupled models, indicating that these models are more slopesensitive than the weakly-coupled model. Stars indicate parameter values used in panels A-D. F: Slope-sensitivity measured with Δg Na , the ratio of threshold g Na values computed from responses to input ramps with slopes 1000 pA/ms and 500 pA/ms. Larger values of Δg Na indicate coupling configurations that are more sensitive to the rising slope of inputs, and are found for models with strong forward coupling and weak backward coupling. https://doi.org/10.1371/journal.pcbi.1006476.g008 Coupling configurations that enhance coincidence detection Comparison of models in the V 2 − h phase plane. We gain an additional perspective on how spiking dynamics of the two-compartment model depend on coupling configuration by examining the V 2 − h phase plane, where h is the sodium inactivation gating variable. For the weak, forward, and strong coupling configurations, we plot h and V 2 -nullclines in Fig 9. The h-nullcline, defined by the function h = h 1 (V 2 ), is shown in black in each panel and is the same for all coupling configurations. Nullclines for V 2 are shown with colored lines in each panel. We calculate these nullclines by setting V 1 to fixed values: -58 mV (V rest ), -40 mV, and -30 mV. These curves are sections (at the selected V 1 values) through the null surfaces of the full three-dimensional V 1 − V 2 − h phase space. The rationale for computing V 2 -nullclines for certain fixed values of V 1 is that V 1 can be viewed as the input to Cpt2 (recall Eq 1). Notice that the V 2 -nullclines are truncated in these figures. The right branches of each V 2 -nullcline, along which the value of h increases with V 2 , are shown in the insets, as are the local maxima (left knee) of the V 2 -nullclines if they exist at values of h between 0 and 1.
We use these phase plane diagrams to illustrate the distinction between tonic firing in the weakly-coupled model and phasic firing in the strongly-coupled and forward-coupled models. For all models, the V 2 -nullclines shift down and to the right as V 1 increases. In the case of weak coupling, the intersection between the h and V 2 -nullclines (the V 1 -dependent fixed point of the system) crosses over from the left branch of the V 2 -nullcline to the middle branch as these nullclines shift. This signals the destabilization of the fixed point via a Hopf bifurcation as V 1 increases beyond a threshold value. For the forward-coupled and strongly-coupled models, by contrast, the left knee (local maximum) of the V 2 -nullcline disappears as V 1 increases. For these models, then, dynamics of the spike-generator cease to be excitable for large V 1 values. This prevents the possibility of repetitive (tonic) firing to steady inputs.
These phase plane diagrams also reveal the role sodium inactivation plays in enhancing slope-sensitivity in models with strong forward coupling. In Fig 9, we observe that increases in V 1 cause larger rightward shifts in the V 2 -nullcline for the strongly-coupled and forward-coupled model, as compared to the weakly-coupled model. As a result, for identical V 1 values, the Coupling configurations that enhance coincidence detection fixed points in the V 2 − h plane are located at larger V 2 values and smaller h values for the models with stronger forward coupling. By construction, there is less attenuation from Cpt1 to Cpt2 in the models with stronger forward coupling: κ 1!2 is larger for the forward-coupled and strongly-coupled models than for the weakly-coupled model. In terms of the passive conductance parameters (recall Fig 3): g 2 is smaller in the forward-coupled model than in the weakly-coupled model, and g c is larger in the strongly-coupled model than the weakly-coupled model. Consequently, a given depolarization of V 1 produces larger depolarizations of V 2 in the forward-and strongly-coupled models compared to the weakly-coupled model. Since larger depolarizations of V 2 cause greater inactivation of sodium current (reducing the h gating variable), inputs that elicit spikes in the strongly-coupled and forward-coupled models must depolarize the cell rapidly to activate sodium current on a timescale faster than inactivation of sodium by the h gating variable. In sum, slope-sensitivity in this model arises from the dynamic, voltage-gated, negative feedback of sodium inactivation, and models with stronger forward coupling are more effective at engaging this process.
Combination of strong forward coupling and weak backward coupling shortens the refractory period. We have discussed phasic firing and slope-sensitivity as factors that make models with strong forward coupling more precise coincidence detectors than models with weak forward coupling. We have not yet, however, established why the particular combination of strong forward coupling and weak backward coupling is advantageous for high-frequency coincidence detection. To do this, we investigate refractory periods and post-spike recovery dynamics. Specifically, we show that models that have both strong forward coupling and weak backward coupling have short refractory periods and rapid recovery after a spike.
The responses in Fig 10 illustrate how refractory period changes with coupling configuration. The stimuli for these simulations are a sequence of two EPSG events (imagine, for example, a brief input from each of the two ears). The amplitude of each EPSG event is three times the size of a unitary synaptic input from the auditory nerve model. The first EPSG in the sequence evokes a spike and, depending on the time delay and the coupling configuration, the trailing EPSG may or may not evoke a spike (even though it is the same amplitude as the first EPSG). In the simulations shown, the weakly-coupled model does not produce a second spike for a time delay of 2 ms, but the forward-coupled and strongly-coupled models do (Fig 10B-10D). In all these simulations, g Na was set g ref Na . To further investigate the effect of coupling configuration on refractory period for these three models, we determined the smallest g Na value for which both EPSGs produced spikes, for varying time delays (Fig 10E). The amplitude of each EPSG in the pair of inputs was three times the amplitude of a unitary auditory nerve fiber input (same as in Fig 10A). This calculation identifies the refractory period-the smallest time delay for which the neuron produces a spike in response to both inputs. Refractory periods in these simulations are 0.7 ms for the forward-coupled model, 0.95 ms for the weakly-coupled model, and 1.1 ms for the strongly-coupled model. In addition to having the smallest refractory period among these three models, the forward-coupled model also has the lowest thresholds across all time-delays. This indicates that the forward-coupled does not require excessive sodium conductance to respond rapidly to successive inputs.
To probe the refractory period across all coupling configurations, we fixed g Na at reference values, and measured the smallest time delay for which the model neuron could respond to both EPSGs. By this measure, models with strong forward coupling and weak backward coupling have the smallest refractory period (lower right corner of Fig 10F).
We explain the effect of coupling configuration on refractory period by examining postspike recovery for the weakly-coupled, forward-coupled, and strongly-coupled models in the V 2 − h phase plane (Fig 11). The monotonically decreasing (black) curve is the h-nullcline and the cubic curves (only partially shown) are the V 2 -nullclines for V 1 fixed at rest and at -42 mV. These are similar to the nullclines shown in Fig 9. The thin (green) curves show trajectories of the response of each model in the V 2 − h phase plane to a pair of synaptic events separated by a 1.5 ms time delay. Sodium conductance is set to g ref Na and the first synaptic event elicits a spike. Responses of the models to the trailing inputs depend on refractory period. For these parameter sets, only the forward-coupled model can respond to both inputs in the pair of successive synaptic inputs.
To gain a qualitative sense of why the forward-coupled model spikes in response to both inputs (and therefore has a shorter refractory period), we compare the state of the neuron as it evolves through the V 2 − h phase plane to the height of the left knee of the V 2 -nullcline. As V 1 increases in response to an excitatory input, the V 2 -nullcline shifts downward (transition from blue to red curves). More specifically, the height of the left knee of the V 2 -nullcline shifts downward as V 1 increases. If an excitatory input (which depolarizes V 1 ) shifts the left knee of the V 2 -nullcline below the position of the trajectory in V 2 − h phase plane, then the V 2 variable will move (quickly, on a fast time-scale) to the right branch of the V 2 -nullcline. This represents the spike upstroke, and is not visible in its entirety in Fig 11 because the right branch of the V 2 -nullcline is outside the field of view of these figures. In these phase plane figures, we see that the second excitatory input, which depolarizes V 1 to roughly -42 mV, only elicits a spike in the forward-coupled model.
There are two factors that contribute to the short refractory period in the the forward coupling model. First, observe that the height of the left knee of the V 2 -nullclines are lower for the models with strong forward coupling (forward-coupled and the strongly-coupled models) than for the weakly-coupled model. In models with strong forward coupling, an increase in V 1 propagates to V 2 with minimal attenuation. This increase in V 2 activates sodium current which lowers the height of the left knee of the V 2 -nullcline. As a consequence, models with strong forward coupling are more excitable (more responsive to an excitatory input) during the recovery period than models with weak forward coupling.
Second, observe that h (the sodium inactivation gating variable) recovers more quickly in models with weak backward coupling. To see this, compare the height of the lower green triangle across the models (lower triangles represent the state of the neuron 1.5 ms after the onset of the first synaptic input). The gating variable has the value h = 0.08 for the weakly-coupled and forward-coupled models at that instant, whereas h = 0.05 for the strongly-coupled model. Sodium inactivation recovers more slowly in the strongly-coupled model. Spike recovery is slow in the strongly-coupled model because the spike in Cpt2 depolarizes V 1 (due to the strong backward coupling), which then prevents V 2 from rapidly returning to rest (due to the strong forward coupling). In models with weak backward coupling, by contrast, V 1 does not depolarize substantially during a spike, and thus Cpt1 can act as a current sink to help return V 2 to rest. These differences in sodium inactivation recovery are subtle, but the combination of the two factors discussed above gives the forward-coupled model a double advantage in responding to high-frequency inputs. Excitation in Cpt1 transfers to Cpt2 efficiently to depolarize the spike-generator, and the spike-generator resets quickly (via recovery of the h gating variable) to accommodate rapid generation of action potentials.

Coincidence detection sensitivity in two-compartment models with dynamic KLT current
In the preceding sections, we have detailed the advantages of strong forward coupling generally, and weak backward coupling for high-frequency stimuli, for coincidence detection sensitivity in a two-compartment neuron model. With the exception of the spike-generating sodium current, the two-compartment model we have considered to this point has been passive. We questioned how our findings would change if additional voltage-gated currents were included. Of particular interest in the context of neural coincidence detection in the MSO is the low threshold potassium (KLT) current. This current is prominent in MSO neurons and enhances their coincidence detection sensitivity [33]. We therefore repeated our test of coincidence detection sensitivity with dynamic KLT conductance (see Materials and methods).
In Fig 12 we show coincidence detection sensitivity measured from responses to three input frequencies (from left to right: 300 Hz, 500 Hz, and 700 Hz). In the top row, 10% of the total conductance in Cpt1 at rest is dynamic KLT conductance. In the bottom row, 10% of the total conductance in Cpt2 at rest is dynamic KLT conductance. The format of each panel is similar to Fig 5 with the color scale in each panel representing the maximal firing rate difference between in-phase and out-of-phase inputs.
Upon comparing these results to our previous results using a passive model (frozen KLT, Fig 5), we observe some differences. Dynamic KLT in the input region (Cpt1, top row) improves coincidence detection sensitivity for all model configurations. While models with strong forward coupling and weak backward coupling (lower right corner of each panel) remain as effective coincidence detectors, the optimal configuration shifts to models with stronger backward coupling. For the 300 Hz stimulus, for instance, the largest firing rate differences are achieved for models with strong forward and strong backward coupling (upper right corner). Models with strong backward coupling can more effectively make use of the KLT current because V 2 spikes propagate back into Cpt1 to activate the KLT current.
Dynamic KLT in the output region (Cpt2, bottom row) also improves coincidence detection sensitivity for all model configurations, but the greatest increases are in models with weak coupling. As a result, coincidence detection sensitivity is nearly uniform across all model configurations, especially in responses to lower and higher frequency inputs. Dynamic KLT in Cpt2 tends to provide the most benefit for models with weak coupling because it provides a secondary source of voltage-gated, dynamic, negative feedback in these models for which sodium inactivation does not suffice to establish dynamics conducive to coincidence detection, including phasic responses to steady inputs and slope-sensitivity (recall Figs 7 and 8).
We compare coincidence detection sensitivity across a range of stimulus frequencies for models with the weak, forward, and strong coupling configurations and that include dynamic KLT, see Fig 13. As above, we test models with dynamic KLT conductance in Cpt1 (Fig 13A1  and 13B1), and models with dynamic KLT conductance in Cpt2 (Fig 13A2 and 13B2). Results for models with dynamic KLT are shown in thick lines. For reference, we also include our earlier results using the frozen KLT model (thin lines, same as results shown in Fig 6). The results are consistent with our observations from Fig 12. In particular, we find that dynamic KLT conductance improves coincidence detection sensitivity (relative to the passive model) for all coupling configurations and nearly all stimulus frequencies. Moreover, the benefit of KLT depends on coupling configuration. Dynamic KLT added to Cpt1 (soma) improves coincidence detection sensitivity the most for the strongly-coupled model. Dynamic KLT added to Cpt2 (axon) improves coincidence detection sensitivity the most for the weakly-coupled model.

Discussion
We systematically examined how soma-to-axon coupling affects neural coincidence detection. We characterized coupling configuration by two parameters (κ 1!2 and κ 2!1 ) representing strength of coupling between input region (soma and dendrite) and output region of a cell (axon and axon initial segment). We also identified a separation of time scales between the slow subthreshold dynamics in the input region and the fast dynamics in the spike generator. We measured coincidence detection sensitivity in the model neurons by viewing them as observers of their own synaptic inputs performing a signal detection task-we interpreted spiking in response to coincident inputs as a measure of correct detection (hits) and spiking responses to non-coincident inputs as a measure of incorrect detection (false positives). Coincidence detection sensitivity in two-compartment models with dynamic KLT. Coincidence detection sensitivity is measured using three different stimulus frequencies (from left to right: 300 Hz, 500 Hz, 700 Hz). Color scheme in each panel is normalized so that the lowest color (dark blue) is the firing rate difference that is 65% of the highest value in that panel. In many panels these values are not reached and range of colors presented does not include these low blue colors. A: 10% of leak conductance in Cpt1 replaced with voltage-gated, low-threshold potassium (KLT) conductance. B: 10% of leak conductance in Cpt2 replaced with voltage-gated, low-threshold potassium (KLT) conductance. https://doi.org/10.1371/journal.pcbi.1006476.g012 Combining these analyses enabled us to elucidate how coupling configuration (described by a few parameters) affects coincidence detection properties of the two-compartment model neuron.
We fixed parameter values based on known properties of principal cells of the medial superior olive (MSO). These neurons are among the first binaural neurons in the mammalian auditory brainstem. MSO neurons receive inputs that originate from both ears and respond Dependence of coincidence detection sensitivity on stimulus frequency for models including dynamic KLT. Firing rate differences for models with dynamic KLT are shown in thick lines. For reference, we also show firing rate differences for models with frozen KLT as thin lines (frozen KLT model, same as Fig 6A). Top row (A1 and B1): 10% of resting conductance in Cpt1 (input region) is dynamic KLT conductance. Bottom row (A2 and B2): 10% of resting conductance in Cpt2 (output region) is dynamic KLT conductance. Left column (A1 and A2): Coincidence detection sensitivity measured as the maximum firing rate difference between responses to in-phase and out-of-phase stimuli. Right column (B1 and B2): Coincidence detection sensitivity measured as the maximum firing rate difference between responses to in-phase stimuli and stimuli with 500 μs time difference. In all panels: coincidence detection sensitivity is measured from responses to a range of frequencies (200 Hz to 700 Hz, in increments of 100 Hz), using the weakly-coupled (blue), forward-coupled (green), and strongly-coupled (red) models. preferentially to coincident inputs. The remarkable temporal sensitivity of these coincidence detector neurons is critical for sound localization processing in mammals.

Strong soma-to-axon coupling enhances neural coincidence detection
Specializations that support temporally-precise coincidence detection in MSO neurons include voltage-gated currents active at membrane potentials near resting values [53], fast and welltimed excitatory synapses [27], and dendritic structure (bipolar dendrites, that segregate inputs from opposite ears onto opposite dendrites) [29,31], see also [25] for review. In this work, we showed that soma-axon coupling is an additional structural specialization that can enhance neural coincidence detection.
By performing a thorough search through the space of coupling configuration, we found that strong forward (soma-to-axon) coupling improved coincidence detection sensitivity. And, moreover, the asymmetric forward-coupled configuration of strong forward coupling and weak backward coupling was the optimal configuration for coincidence detection in response to higher frequency inputs (500 Hz to 700 Hz). We identified advantages of strong forward coupling for neural coincidence detection including phasic and slope-sensitive spiking dynamics, and (for the forward-coupled model) short refractory periods. These advantages depended on the action of the sodium inactivation gating variable, which is the sole source of voltage-gated, dynamic, negative feedback in the version of the two-compartment model with frozen KLT.

Dynamic KLT current further enhances neural coincidence detection
KLT current is an additional source of negative feedback and one known to be prominent in MSO neurons in soma and dendrites regions [31], as well as axon regions [38]. We found that dynamic KLT current improved coincidence detection for nearly all coupling configuration and stimulus frequencies. There were notable interactions between coupling configuration and KLT current. Dynamic KLT current improved coincidence detection in neurons with strong forward and backward coupling so that this strongly-coupled configuration (the configuration most similar to a one-compartment model) became optimal for coincidence detection in response to intermediate-frequency stimuli. In addition, dynamic KLT current localized to Cpt2 (the spike-generator region of the neuron model) could rescue coincidence detection sensitivity in neurons with weak soma-axon coupling so that this weakly-coupled configuration could exhibit coincidence detection sensitivity on par with other models. There was some loss of efficiency when dynamic KLT current was included (g ref Na values increased by 5% to 30%), but these differences are modest compared to the order of magnitude differences in g ref Na across coupling configurations.
Our work adds to the characterization of MSO neurons as coincidence detector specialists equipped with an array of features-structural and dynamical-that enhance their temporal precision. Strong forward coupling appears to be a natural configuration for coincidence detection. Moreover, the specific combination of strong forward coupling and weak backward coupling is advantageous for high frequency coincidence detection. The benefits of these coupling configurations can be supplemented with appropriately-targeted KLT current. The need for multiple, complementary mechanisms that enhance coincidence detection MSO neurons has been explored previously, but usually for one-compartment models (Na inactivation and KLT activation as two sources of negative feedback [36]; KLT current and hyperpolarizationactivated cationic current as currents that regulate input resistance [32]). The effect of KLT is considerable and well-studied, so here we emphasized the role of structure (soma-to-axon coupling). The intrinsic advantage of the forward-coupling configuration may help maintain coincidence detection sensitivity in scenarios in which KLT is less effective(early in development [73], for instance, or when KLT is inactivated by long time-scale channel gating dynamics that are not in our model [32]).
Other features of MSO neurons such as inhibitory synaptic inputs [57,[59][60][61][62][63] and dendritic structure [29][30][31] specialize these cells for coincidence detection. The variety of physiological tools MSO neurons use to perform coincidence detection emphasizes the exceptional nature of the temporally-precise computations these neurons perform. Inhibitory inputs in MSO primarily target cell bodies [74] so careful study of the inhibition, excitation, and cell-structure would be an interesting avenue of future research. Previous studies using point-neuron models and dynamic clamp electrophysiology demonstrate that the interaction of inhibition and KLT current can shift the location of the peak of time difference tuning curves [57,61,63]. We would expect similar effects in a two-compartment model, as long as the soma-targeting inhibition and KLT current are not electrically isolated. These effects may require, therefore, a prominent KLT current in the soma (co-localized with inhibition) or a strong soma-to-axon coupling configuration to enable soma-targeting inhibition to interact with axonal KLT current (if it is present).
We tested coincidence detection sensitivity with out-of-phase inputs and inputs with a fixed 500 μs delay. We observed qualitatively similar results for both types of inputs. The latter stimulus may be more relevant in studying neural coincidence in the context of sound localization. Time differences in this context would be created by differences in travel times of sounds arriving at the two ears (interaural time differences) and are limited by animal head size. In humans, for instance, maximal interaural time differences created by head size are approximately 700 μs.

MSO neurons appear forward-coupled
Our model is phenomenological-the two-compartments are lumped representations of input and output regions. We do not, therefore, resolve structural details of dendrites or spike initiation zones (see [38] for an example of the latter). Nonetheless, we can make qualitative observations that relate our findings to MSO physiology.
Action potentials in the MSO are likely generated in a spike initiation zone near the soma, and back-propagated action potentials in the soma are small and graded [37]. This indicates a strict electrical segregation of the soma and dendrites from the axonal initiation zone, (in the words of [37]). In the context of our model, this corresponds to weak backward coupling (small value of κ 2!1 ). Backpropagation of signals into the dendrites is further attenuated due to the low input resistance of these neurons and the strong effects of voltage-gated potassium current [31]. Additionally, current injection into the soma reliably evokes action potentials that propagate into the axon [37]. This suggests a configuration in which the soma has a strong effect on the spike-generator (minimal attenuation, large value of κ 1!2 in our model).
Taken together then, it appears that MSO neurons may be structured in a forward-coupled manner, consistent with our observations that this configuration confers advantages for coincidence detection by engaging sodium inactivation as dynamic negative-feedback mechanism, by promoting rapid resetting of the spike generator (shortening the refractory period, which enables high-frequency spiking), and by enabling efficient spike generation (smaller sodium conductance required).
The complete picture of MSO excitability and axonal structure is, of course, more complicated than our two-compartment description. Recent computational simulations provide evidence that spike generation may occur throughout MSO axons (initial segment and multiple nodes of Ranvier) [38]. Spike generation at more distal locations on the axon can preserve excitability in response to high-frequency stimuli by preventing inactivation of sodium channels [38,44]. Studies of coincidence detector neurons in related structures in the avian auditory show that excitability of these neurons can be adjusted via modulation of ion channel density in spike generator regions [42][43][44]. This raises intriguing questions about plasticity in the spike initiation zone, and dynamic regulation of the soma-to-axon connection.

A framework for investigating neural structure, dynamics, and computation
We have formulated a family of two-compartment models to investigate neural coincidence detection in MSO neurons. We showed that parameters in this two-compartment framework can be chosen in a principled manner to explore the range of coupling configurations, while maintaining similar passive dynamics in the input region. With this approach, we identified how structure (the nature of soma-axon coupling) affected dynamics in the spike-generator region, and, in turn, how these differences in dynamics affect the sensitivity of coincidence detector neurons to synaptic inputs.
Our approach provides a unifying view of structure and function in neurons performing an identified computation. It is one that should find applications in studies of other neurons. Coincidence detector neurons in the auditory brainstem of owls, for instance, have been modeled as two-compartment structures [41]. The two-compartment idealization has also been useful for investigating dynamics of bursting [45,49,75], bistability [46], oscillations [47], and resonance [50] in neurons, and could also describe signaling between a (large) dendrite region and a (small) dendritic spine. Our framework for creating, and systematically exploring, a parameter space of soma-axon coupling configurations can be used to shed further light on the relationship between structure, dynamics, and function in these and other neural systems.
Supporting information S1 Fig. Simulations of ROC curves and identification of best g Na . A: Firing rate to in-phase inputs (y-axis) plotted against firing rate to out-of-phase inputs (x-axis) with g Na varied to explore a range of firing rates. B: Firing rate differences (y-axis) plotted against Na conductance (g Na normalized to reference g Na values, see text for details). We identify the value of g Na that is best for coincidence detection by finding the g Na value that results in the greatest firing rate difference. Best g Na and the corresponding firing rate values are marked with dots in these figures. Error bars are standard error of mean firing rates computed from 100 repetitions of 1-second long simulations for each g Na and frequency. We show results for the weakly-coupled, forward-coupled, and strongly-coupled models (colored lines) and responses to three stimulus frequencies (from left to right: 300, 500, and 700 Hz). (TIF) S2 Fig. Reference g Na for models with 10% dynamic low-threshold potassium (KLT) conductance. A: Reference g Na across parameter space of coupling strengths, for model without dynamic KLT. This panel reproduces Fig 2B, see text for definition of reference g Na . B: Reference g Na for two-compartment models with 10% of leak conductance in Cpt1 replaced by dynamic KLT conductance. C: Reference g Na for two-compartment models with 10% of leak conductance in Cpt2 replaced by dynamic KLT conductance. (TIF)

S3 Fig. Effect of stimulus frequency and coupling configuration on best Na conductance.
The value of g Na at which the two-compartment model achieves its best coincidence detection sensitivity (maximal firing rate difference between response to in-phase and out-of-phase inputs) for (A) 300 Hz stimuli, (B) 500 Hz stimuli, and (C) 700 Hz stimuli. D: Detailed view of these best g Na values for the weakly-coupled, forward-coupled, and strongly-coupled models as function of input frequency. (TIF) S4 Fig. Coincidence detection tuning curves for the weak, forward, and strong coupling models. Firing rates in response to 300 Hz (left column), 500 Hz (middle column) and 700 Hz (right column) stimuli. Sodium conductance values are selected so that the firing rate for coincident inputs (0 ms time difference) for each model is 150 spikes/second (top row), 250 spikes/ second (middle row), or 350 spikes/second (bottom row). (TIF) S1 File. Excel spreadsheet of figure data. (XLSX)