Abstract
Extracellular field potentials (EFPs) are an important source of information in neuroscience, but their physiological basis is in many cases still a matter of debate. Axonal sources are typically discounted in modeling and data analysis because their contributions are assumed to be negligible. Here, we established experimentally and theoretically that contributions of axons to EFPs can be significant. Modeling action potentials propagating along axons, we showed that EFPs were prominent in the presence of terminal zones where axons branch and terminate in close succession, as found in many brain regions. Our models predicted a dipolar far field and a polarity reversal at the center of the terminal zone. We confirmed these predictions using EFPs from the barn owl auditory brainstem where we recorded in nucleus laminaris using a multielectrode array. These results demonstrate that axonal terminal zones can produce EFPs with considerable amplitude and spatial reach.
https://doi.org/10.7554/eLife.26106.001Introduction
Extracellular field potentials (EFPs) are at the heart of many experimental approaches used to examine the inner workings of the brain. Types of EFPs include invasively recorded signals such as the electrocorticogram (ECoG) and the local field potential (LFP), as well as the noninvasively recorded electroencephalogram (EEG) and the auditory brainstem response (ABR) (Brette and Destexhe, 2012; Nunez and Srinivasan, 2006). Measures derived from the EFP such as the current source density (CSD) and multiunit activity (MUA) are also frequently used. The origins of these signals and measures, especially in cases in which the activity is not clearly attributable to a single cell, is a matter of debate (Buzsáki et al., 2012).
EFPs in the brain were long thought to be primarily of synaptic origin (Buzsáki et al., 2012). As a consequence, many modeling studies focussed on the extracellular fields induced by postsynaptic currents on the dendrites and the soma of a neuron (Holt and Koch, 1999; Einevoll et al., 2013; Lindén et al., 2011; Lindén et al., 2010; FernándezRuiz et al., 2013). However, a number of recent data analyses and modeling efforts have revealed that active, nonsynaptic membrane currents can play an important role in generating populationlevel EFPs (Reimann et al., 2013; Anastassiou et al., 2015; Schomburg et al., 2012; Ray and Maunsell, 2011; Belluscio et al., 2012; Waldert et al., 2013; Ness et al., 2016; SchefferTeixeira et al., 2013; Reichinnek et al., 2010; Sinha and Narayanan, 2015; Taxidis et al., 2015), including far reaching potentials detectable at the scalp (Baker et al., 2011; Teleńczuk et al., 2015). Currents from the axon are still thought to be so small as to be of minor importance for the EFP.
One of the reasons for the assumption that axonal currents contribute little to the EFP is that the far field of an action potential traveling along an idealized straight and long axon is quadrupolar, meaning that it decays faster with distance than synaptic sources, which are typically dipolar (Nunez and Srinivasan, 2006). Surprisingly, theoretical (Tenke et al., 1993) and experimental studies indicated that the EFP of axonal responses may also have a dipolar structure. For example, Blot and Barbour (2014) reported an EFP with a characteristic dipolar structure in the vicinity of cerebellar Purkinje cell axons; other studies (Swadlow et al., 2002; Swadlow and Gusev, 2000) showed that the axonal part of the EFP of thalamocortical afferents showed a polarity reversal associated with a dipolar source, and classical current source density studies found dipolar current distributions in axonal terminal zones in the cortex and the lateral geniculate nucleus, and attributed these to axons because of conduction velocities (Mitzdorf and Singer, 1978; Mitzdorf, 1985; Mitzdorf and Singer, 1977; see also Tenke et al., 1993). Here we introduce another experimental system and show a strong dipolar, axonal field potential in the auditory brainstem of the barn owl.
The discrepancy between the quadrupolar structure of EFPs generated by idealized axons, and the experimentally observed dipolar structure raises the question of how axons are able to generate dipolar field potentials. In this manuscript we show how dipolar far fields in the EFP of axons can be explained by the axons’ anatomical structure. In particular, the branchings and terminations of axons in their terminal zone area deform the extracellular waveform (Gydikov et al., 1986; Gydikov and Trayanova, 1986; Plonsey, 1977) and can lead to a dipolar EFP structure. Axon bundles, sometimes called fascicles, exist throughout the peripheral and central nervous system and have such terminal zones (Kandel et al., 2000; Hentschel and van Ooyen, 1999; Nornes and Das, 1972; Goodman et al., 1984). The white matter of the mammalian brain can be viewed as an agglomeration of such fascicles (Schüz and Braitenberg, 2002). We therefore predict pronounced contributions of axon bundles to EFPs, which are neglected in current models.
In what follows, axonal contributions to the EFP are first investigated by a numerical model based on forward simulation (Holt and Koch, 1999; Gold et al., 2006). This first model includes a largescale multicompartment simulation (Jack et al., 1975; Rall, 1959; Abbott, 1992; Hines and Carnevale, 1997; Hines et al., 2009) of an axon population. We then outline the basic mechanisms by means of a second, analytically tractable, model of a generic axon bundle. Finally, we validate model predictions with data from multisite invivo electrophysiological recordings from the barn owl auditory brainstem.
Results
Effects of axonal bifurcations and terminations on extracellular action potentials
To understand how the geometry of an axon affects the extracellular waveform associated with action potentials, we first numerically simulated single action potentials propagating along generic axons and calculated their contribution to the EFP (for details, see Materials and methods). This was done for five scenarios: quasiinfinite axons, terminating axons, bifurcating axons, axons that bifurcate as well as terminate, and axon bundles (Figure 1). We began by simulating a long axon, approximating an infinitely long axon following a straight line path, neither bifurcating nor terminating (Figure 1A). The extracellular action potential has the characteristic triphasic shape. As the action potential travels along the axon, the waveform is translated in time with the conduction velocity, but is otherwise unchanged. The triphasic shape is also present in the spatial arrangement of transmembrane currents at any given time, which is the reason for the quadrupolar EFP response traditionally assumed for axons.
There are two ways of understanding the triphasic shape of the extracellular waveform. One way is by attribution of the peaks of the response to individual current types. The first, small positive peak corresponds to the capacitive current, the large and negative second peak to the sodium current, and the final positive peak to the potassium current (Gold et al., 2006). Another, more mathematical way of understanding the triphasic shape is specific to the nature of the axon. Due to Kirchhoff’s law and cable theory (see Materials and methods), the local transmembrane current in a homogeneous axon is proportional to the second spatial derivative of the membrane potential along the direction of the axon. Because the action potential is roughly a traveling wave, the currents are also proportional to the second temporal derivative of the membrane potential. The three extrema of the EFP are thus related to the points of maximum curvature in the action potential waveform, namely the onset, the maximum, and the end of the spike.
Next we simulated the response of an axon that terminates (Figure 1B). Here the action potential approaching the recording location (top traces) has the same, triphasic EFP response as in the nonbranching case. When the action potential reaches the termination point, its EFP gradually deforms into a biphasic response, with a positive peak preceding a negative peak. The mechanism for this deformation can be understood as follows: As the action potential approaches the recording location next to the termination, the majority of the transmembrane currents flow at points located before the termination, and they are almost identical to those in the nonterminating case; the first, capacitive peak is not affected. As already mentioned, the second and third peaks of the extracellular action potential in the nonterminating case are generated by currents close to or after the electrode location. In the terminating axon, there are no currents at points after the termination, leading to a partial suppression of the second peak and a complete suppression of the third peak.
Another generic structure found in axons is a bifurcation. To emphasize the impact of bifurcations, we simulated a single axon that bifurcates three times on each branch within a distance of 200 µm (100 µm between branchings), leading to a total number of 8 collaterals leaving the bifurcation zone (Figure 1C). (Note that in order to avoid confounding effects, the horizontal distances between axons in Figure 1C–E are for illustration only; all collaterals were simulated to lie on a straight line.) The EFP far away from the bifurcation zone has a triphasic shape and resembles the one observed in Figure 1A, and the amplitude is proportional to the number of axon fibers. The EFP near the bifurcation zone has a biphasic shape. Although there is an initial tiny positive peak, the response is dominated by the second, negative and the third, positive peak. This waveform can again be understood by comparison to the first example (Figure 1A) containing the infinitely long axon: The tiny positive initial peak resembles the infinite case, because it is constituted by the action potentialrelated currents flowing within the part of the axon before the bifurcation. As the action potential passes the bifurcation zone, there are now several action potentials (one in each fiber). Because of the active nature of the action potential, the active currents are the same in each outgoing fiber as in the incoming fiber. This leads the second and third peak to be multiplied in size, yielding a quasibiphasic response. We chose to simulate several bifurcations because this leads to a clearer effect in the EFP. In the case of a single bifurcation, this effect is also present, but the amplification of the second and third peak relative to the first peak is not as notable as in this example.
To understand how bifurcations and terminations interact when they are present in the same axon, we simulated an axon with an identical number of bifurcations as in the previous case, but then added terminations to all the fibers 700 µm after the bifurcation zone (Figure 1D). We found that this configuration leads to the same biphasic responses as observed in the cases of the isolated anatomical features. A triphasic response occurred inbetween the bifurcation and termination zones. A notable point here is that the potential at the bifurcation and termination are both biphasic and on the same timescale, but opposite in polarity.
After having studied the EFPs of single axons, we next started to simulate axon bundles, because axons often run in parallel bundles in the brain. Moving towards more biologically plausible axon geometries, we considered bundles consisting of axons with slightly altered spatial arrangement: we randomly perturbed the precise locations of bifurcations and terminations in the axon tree (Figure 1E, left; for details, see Materials and methods). We simulated 100 axons and stimulated each axon with an inhomogeneous Poisson spike train (Kuokkanen et al., 2010; Softky and Koch, 1993). The driving rate of the inhomogeneous Poisson process was the same for all axons and consisted of a constant background rate (100 spikes/s) and a Gaussian pulse of heightened activity (2000 spikes/s). The standard deviation of the pulse was 1 ms, resulting in an additional 3.5 spikes per axon being fired on average over the entire duration of the pulse. The resulting extracellular populationlevel waveforms (Figure 1E, right) showed a polarity reversal reminiscent of Figure 1D. However, in the bifurcation zone, the summed contribution from many fibers and action potentials lead to a monophasic negative peak, and in the termination zone there was a monophasic positive peak. Interestingly, the summed potential at the center of the terminal zone largely cancelled out.
The fact that the responses in Figure 1E were mostly monophasic can be explained by the presence of a nonzero bias in the biphasic responses observed for the single spike responses in Figure 1D: close to a bifurcation, the area under the negative part of the curve slightly exceeded that of the positive part, and vice versa close to a termination. When summed up over many spikes with different timings, this difference in areas induced a positive or negative polarity of the population response in Figure 1E.
The reversal behaviour shown in Figure 1E is similar to the polarity reversal associated with a dipole observed in experimental studies (Swadlow and Gusev, 2000; Swadlow et al., 2002; Blot and Barbour, 2014). To summarize, simple onedimensional model axon structures can produce complex and diverse spatiotemporal EFP responses, including monophasic, biphasic and triphasic waveforms, comparable to experimentally recorded responses.
Axonal projections generate a dipolelike field potential
Dipolelike EFPs have a much larger spatial reach than quadrupolarlike EFPs, which are typically associated with axons (Nunez and Srinivasan, 2006). To further understand whether and how axons can generate a dipolar EFP, in Figure 2 we turned to threedimensional axon morphologies, in contrast to the onedimensional case studied in Figure 1 (for details, see Materials and methods). We thus simulated a parallel fiber bundle of 5000 axons that at first runs at a constant number of fibers without bifurcations and then reaches a terminal zone. Within this terminal zone, the fibers first bifurcate, which increases the number of fibers. Finally, as the axons reach the end of the terminal zone, they terminate and the number of fibers decreases to zero (example axon shown in Figure 2A). To model more closely the actual axonal structures found in nature, we included a radial fanning out of the branches as well as a more diverse set of morphologies with a variable number of bifurcations and terminations (for details, see Materials and methods).
The spiking activity of a generic axon bundle was simulated by a background spontaneous firing rate of 100 spikes/s and a short pulse of increased activity. We chose a Gaussian pulse with an maximum instantaneous rate of 2900 spikes/s and a standard deviation of 2.8 ms. Note that this high driving rate is only the instantaneous maximum, and the actual firing rate is limited by the refractory period following a spike. These numbers are motivated by the early auditory system of barn owls (Köppl, 1997a; Sullivan and Konishi, 1984; Konishi et al., 1985), where instantaneous spike rates of 3000 spikes/s occur in response to click stimuli (Carr et al., 2016). However, our approach is not limited to the auditory system (which would also require the introduction of the synchronization of the spike times to the auditory stimulation frequency, called phase locking). Instead, this pulse of activity could relate to various kinds of evoked activity in the nervous system, such as sensory stimulation, motor activity or a spontaneous transient increase in population spiking activity.
To characterize the spatiotemporal dynamics of the evoked EFP, the time course of the potential was calculated for several locations along the axon trunk. The responses were averaged over 10 repetitions. We divided our analysis into two frequency bands by filtering the responses. The first frequency band was obtained by lowpass filtering with a cutoff frequency of 1 kHz (Figure 2A–C). The second frequency band was the multiunit activity (MUA) obtained by highpass filtering with a cutoff frequency of 2.5 kHz (Figure 2D–F). To make the MUA easier to interpret in terms of overall activity reflected, it was halfwave rectified and lowpass filtered (<500 Hz, see Materials and methods). The two frequency bands showed a qualitatively different spatiotemporal response in the vicinity of the projection zone, as we will show in the following.
We first studied the effect of the Gaussian rate pulse on the lowpass filtered EFP (Figure 2B). The filtering removed most of the identifiable components of individual spikes, while a populationlevel signal remained, similar to an LFP signal. The distribution of the maximum amplitudes of these responses is shown by the colored contour lines in Figure 2A and the colored voltage traces in Figure 2B. Surrounding the terminal zone of the axon bundle in Figure 2A, lowpass filtered EFP amplitudes showed a doublelobed shape typical of a dipole.
In Figure 2B, the lowpass filtered EFP responses mostly showed monophasic deflections elicited by the population firing rate pulse, in a manner similar to that observed in Figure 1E. Such deflections were visible at all recording locations. In the radial direction away from the axon tree, i.e. in the horizontal direction in the figure, the lowpass filtered EFP amplitude decays. In the axial direction along the axon tree, i.e. in vertical direction in the figure, the voltage deflection reverses polarity in the middle of the terminal zone of the bundle (Figure 2B). The polarity reversal occurs by a decrease of the amplitude to zero and a subsequent reappearance with reversed polarity (as opposed to a polarity reversal through a gradual shift in phase). This behaviour is also typical for a dipolar field potential.
The point of the polarity reversal coincides with the middle of the terminal zone. Interestingly, this means that the absolute value of the response amplitude reaches a local minimum just at the axial location at which the number of axonal fibers reaches a maximum. To better understand how the anatomical features of the axon bundle and the EFP response amplitude are related, we compared its signed maximum value (meaning the signed value corresponding to the maximum magnitude of the amplitude of the EFP, black line in Figure 2C) with the change of the number of nodes per 200 µm bin (purple histogram in Figure 2C): Along the nerve trunk the number of fibers is constant. As the axon bundle reaches its terminal zone, the number of bifurcations increases (purple bars point to the right in Figure 2C). The increase of bifurcations is followed by an increase in terminations. In the middle of the terminal zone, the number of bifurcations and terminations are equal. At the same depth, the amplitude of the EFP component crosses zero. At the end of the terminal zone, the terminations outweigh the bifurcations (purple bars point leftwards in Figure 2C). As the axon bundle ends, there are no longer any bifurcations or terminations, and the number of fibers decays toward zero. Overall, the signed maximum amplitude EFP (black trace) follows the distribution of branchings and terminations (purple histogram). This progression of amplitudes in the lowfrequency components seen in Figure 2C is also visible in Figure 2B, most clearly in the first column. The polarity reversal in the center of Figure 2B corresponds to the crossing of zero amplitude in Figure 2C.
To understand how the EFP contributions are related to individual spikes, we next turned our attention to the highfrequency MUA. The MUA is thought to reflect local spiking activity (Stark and Abeles, 2007). In Figure 2D, the isoamplitude lines of the MUA appeared like an ellipsoid centered at the terminal zone (Figure 2D); they did not show the doublelobe shape observed for the lowpass filtered EFP in Figure 2A.
The shape of the MUA response was weakly dependent on the recording location. The main change across locations was in the scaling of the amplitude (Figure 2E). The amplitude decays with radial distance from the trunk. In the axial direction, the amplitude reaches its maximum in the middle of the fiber bundle. This dependence of the MUA amplitude on the axial location is further examined in Figure 2F. The amplitude of the MUA component (black trace) changes in accordance with the local number of nodes per unit length (tealcolored histogram), which is proportional to the number of fibers. The local number of fibers and the MUA amplitude are both constant along the nerve trunk. Both measures then increase in amplitude as the number of fibers is increased by bifurcations. As the fibers terminate and the number of fibers decreases, so does the amplitude of the MUA.
To conclude, we have shown a qualitatively different behaviour in the low and highfrequency components of the EFP, i.e. for the lowpass filtered EFP and the MUA. The particular branching and terminating structure of the axon bundle may thus give rise to a dipolar lowpass filtered EFP.
Effects of bifurcations and terminations on distance scaling of EFPs
To further demonstrate that bifurcations and terminations of axons give rise to a dipolar field, we investigated the effect of an axon terminal structure on the spatial reach of the EFP (Figure 3). Motivated by the fundamentally different spatial distributions of the lowpass filtered EFP and the highfrequency MUA in Figure 2, we again differentiated between these frequency bands and simulated an axon bundle containing a terminal zone with bifurcations and terminations. Moreover, as a control, we also simulated an axon bundle without bifurcations in which a fixed number of fibers simply terminates.
In order to separate the effects of any radial fanning out of the axon bundle from the effects of bifurcations and terminations, and to afford better analytic tractability, we transitioned back to a simpler onedimensional model of the axon bundle (see Materials and methods). This model omitted the radial fanning out of the bundle in the terminal zone, as in Figure 1. Furthermore, we discarded the detailed conductancebased simulation of the membrane potential, and instead assumed a fixed membrane potential waveform traveling along the axon trunk with a constant propagation velocity. Using linear cable theory, it was then possible to calculate the membrane currents necessary for the determination of the EFP. The analytic nature of the simplified model also allowed us to consider a continuous number of fibers instead of simulating discrete bifurcations and terminations. All following simulations are based on this simplified model.
To verify that this simpler onedimensional model accurately captures the EFP response of an axon bundle, we applied parameters equivalent to those used in Figure 2 and compared the resulting EFP to that obtained from the full biophysical model. We calculated the relative difference of the EFPs by taking the absolute value of their differences, and normalizing by the sum of their absolute values. Averaged over time, this relative difference at distances greater than 1 mm from the center of the projection zone was <0.05 in axial direction. In the radial direction, we found, as expected, larger relative discrepancies of <0.3 for radial distances > 1 mm. In what follows, we focus on the axial direction, which is the dipole axis.
Let us now specify how we simulated the two axon bundle morphologies. The control case was a nonbifurcating bundle, which had a constant number of 50 fibers up to the termination point, and then tapered out with a Gaussian profile that was centered at the termination point with a height of 50 fibers and width (standard deviation) of 300 µm. The second case was that of an axon bundle with a projection zone containing bifurcations. Here we added to the distribution of the number of axons used for the nonbifurcating control case a further Gaussian distribution to account for the projection zone. This additional Gaussian was also centered at the termination point, but had an amplitude of 450 fibers and a standard deviation of 500 µm, meaning that the overall width of the terminal zone was $\approx $2 mm. Unlike the taperingout in the control condition, this component was added both before and after the termination point. It resulted in a maximal fiber number of 500 at the termination point, which is a factor 10 larger than in the control case. Both distributions constructed in this way were smooth, and they had smooth first derivatives in space. In both cases, the number of fibers decreased monotonically after the termination point. We considered a conduction velocity of 1 m/s in this example, though results are qualitatively the same for other values. In order to understand the frequencyspecific effects of the projection zone, we calculated the responses to membrane potential components with temporal frequencies between 25 Hz and 5 kHz, with the same amplitude for each frequency component. For the conduction velocity 1 m/s, these temporal frequencies corresponded to spatial wavelengths from 10 mm to 0.2 mm, i.e. from much larger to much smaller than the width of the terminal zone. We then calculated for each frequency/wavelength the average amplitude of the resulting EFP response by taking its standard deviation. Due to the linear nature of our model, the frequency responses obtained in this manner are applicable to the Fourier components of any membrane potential timecourse.
The dipolelike component observed in Figure 2 for the lowfrequency component had its dipole axis aligned with the axon trunk. We therefore considered the distance $r$ beyond the termination point in the direction extending the axon trunk, which we called the axial direction. Because we suspected a dipolar response, we expected the amplitude of the field potential to decay as ${r}^{2}$. To test the scaling behaviour of this component, we first plotted the average amplitude of the lowfrequency responses ($f$<1 kHz) in axial direction in Figure 3A. The plot is on a double logarithmic scale, meaning that the slope of the curve corresponds to the scaling exponent, and the vertical offset corresponds to the amplitude of the size of the dipole moment, which is a measure for the strength of the dipolar EFP. We observed the expected ${r}^{2}$ scaling for distances larger than the extent of the bifurcation zone ($\gtrsim 1$ mm).
Comparing the responses of the bifurcating axon bundle and the nonbifurcating control condition (full and dashed lines in Figure 3A), we saw that for short distances (<1 mm) the response of the bifurcating case was a factor 10 larger than the control. At these distances the response was due to the local fibers, of which there are 10 times more in the bifurcating case. At distances larger than 1 mm, we observed that the distance scaling was proportional to ${r}^{2}$, meaning that there was a dipole moment in both conditions (a vanishing dipole moment would have implied a slope steeper than ${r}^{2}$). Interestingly, for distances larger than 1 mm the response in the bifurcating case exceeded the control by a factor 20. We thus concluded that at low frequencies, the bifurcation zone contributes supralinearly to the dipole moment.
The reason for this supralinearity is that contributions from different parts of the axon bundle can interfere constructively or destructively. The maximum constructive interference occurs when the spatial width of the oscillation agrees with that of the projection zone (Figure 3A, top). Importantly, currents from fibers inside the projection zone on average interact destructively with those from outside the projection zone. The magnitude of this effect depends on the ratio of the number of fibers inside the projection zone compared to the number outside. The larger the ratio the smaller the impact of destructive interference. Thus, bifurcations suppress the destructive interference (Figure 3A, bottom, full line). On the other hand, for a ratio of one, i.e. the nonbifurcating case, destructive interference is strong, which diminishes the overall response amplitude (Figure 3A, bottom, dashed line).
Next, we examined the highfrequency (>2.5 kHz) component (Figure 3B). As in the lowfrequency case, the response at distances < 1 mm was greater in the bifurcating case by a factor of 10. The asymptotic scaling was also ${r}^{2}$ for axial distances $>1$ mm in both cases. However, unlike in the lowfrequency case, the amplitudes were similar between bifurcating and nonbifurcating cases. Thus, the presence of a bifurcation zone did not contribute to the highfrequency dipole moments in the EFP. This feature is explained by the small spatial wavelength of the stimulus compared to the width of the bifurcation zone (Figure 3B, top)
Frequencydependence of dipolar axonal EFPs
We showed that the dipole moment depends on both the anatomy, that is, the presence of a projection zone, and the temporal frequency range (low vs. high frequencies) of the underlying activity. This relationship can be qualitatively understood by considering that in an axon bundle a voltage waveform propagates at some conduction velocity. The temporal frequency of this signal thus corresponds to a spatial frequency. If the spatial frequency of the membrane potential matches the width of the projection zone, the dipole moment can reach its maximum. In this case, at some point in time, positive membrane currents flow in one half of the projection zone and negative membrane currents flow in the other half (Figure 3A, top). For example, if the voltage waveform has a temporal frequency of 1 kHz and propagates at 1 m/s, the spatial wavelength is 1 mm. If the spatial width of the termination zone is about 1 mm, the dipole moment is maximal. In contrast, if the spatial wavelength is much smaller than the width of the projection zone, an alignment between projection zone and current flow is not possible, and the dipole is not amplified (Figure 3B, top) (for a detailed derivation see Materials and methods).
To quantitatively understand the frequencyspecific contributions to the dipole moments, we examined the scaling behaviour of the EFP as a function of frequency. The amplitude of the dipole moment was determined by fitting a straight line with slope $2$ to the double logarithmic scaling of the standard deviation of the EFP at a given frequency. The fit was performed for distances $>1$ mm. The extrapolation of this straight line to the axial distance 1 µm was then proportional to the dipole moment.
The normalized frequencyspecific dipole moments are shown in Figure 3C. The dipole moment of the bifurcating case (solid line) has a maximum at around 500 Hz, as expected due to the agreement of the spatial wavelength (1 m/s / 500 Hz = 2 mm; gray vertical line in Figure 3D and C) and axial width (about 2 mm; gray box in Figure 3A and B) of the projection zone. The match is not exact because the shapes of sine wave and Gaussian are different. For lower and higher frequencies, there is a mismatch in spatial wavelength and the width of the projection zone, meaning that the projection zone contributes less to the dipole moment, as also observed in Figure 3B for higher frequencies. In the nonbifurcating control case (dashed line) there is no projection zone, and the dipole moment decays monotonically with rising frequency because higher frequencies correspond to a smaller spatial separation of positive and negative currents, and thus to a smaller dipole moment.
In the bifurcating case, the maximum number of fibers was increased by a factor of 10. Accordingly, an increase in the dipole moment by a factor of 10 from the nonbifurcating to the bifurcating case could be explained by just linearly summing the dipole moments of individual fibers. An increase in the dipole moment by a factor greater than 10 would be supralinear. In Figure 3D we compared this relative impact of the terminal zone on dipole moments (red line) by plotting the dipole moment ratios across frequencies. The contribution of the terminal zone is greater than 10 (dotted line) for intermediate frequencies between about 200 and 1300 Hz, and smaller than a factor 10 outside this frequency range.
Together, these observations show us that the terminal zone makes a frequency specific contribution to the farreaching dipole field potential of the axon bundle. This provides a deeper understanding of the findings of Figure 2: At low frequencies ($<1$ kHz), we observed a supralinear dipolar behaviour due to the specific morphology the bundle, with the projection zone forming the dipole axis. At higher frequencies, the bifurcation zone does not amplify the dipole moment, meaning that we could observe responses mainly locally.
The barn owl neurophonic potential in nucleus laminaris as an example for a dipolar field in an axonal terminal zone
To test our prediction of dipolar extracellular field potential responses due to axon bundles, we recorded EFP responses from the barn owl auditory brainstem. The barn owl has a highly developed auditory system with a strong frequencyfollowing response in the EFP (up to 9 kHz, [Köppl, 1997b]), called the neurophonic, which can be recorded in the nucleus laminaris (NL). In NL, the input from the two ears is first integrated to calculate the azimuthal location of a sound source, and this information is encoded in the EFP (Carr and Konishi, 1990). The EFP in this region is mainly due to the afferent activity, and the contribution of postsynaptic NL spikes is small (Kuokkanen et al., 2010; Kuokkanen et al., 2013). Furthermore, the anatomy of the afferent axons is well known and follows a stereotypical pattern (Carr and Konishi, 1988; Carr and Konishi, 1990): Two fiber bundles enter the nucleus, with fibers from the contralateral ear entering ventrally, and from the ipsilateral ear entering dorsally. The axon bundles reach the NL from their origin without bifurcating, then bifurcate multiple times at the border of the NL, and then terminate within NL. Axon bundles have a strong directional preference and run roughly in parallel. Most of the volume within NL consists of incoming axons. This well studied physiology and anatomy makes the system an ideal candidate to investigate the EFPs of axon bundles; see the Discussion for arguments why synaptic contributions to the EFP could also be neglected here.
To explore the spatiotemporal structure of the EFP in NL, we performed simultaneous multielectrode recordings of the response in NL (Figure 4A) to contralateral monaural click stimuli. The click responses showed distinct lowfrequency (Figure 4B) and highfrequency (Figure 4C) components, as previously reported (Wagner et al., 2009). The frequency of the highfrequency ringing corresponds to the recording location on the frequency map within NL, and the ringing reflects the frequency tuning and phase locking of the incoming axons. In addition, there is a lowfrequency component in the response (Figure 4B). We filtered the data to separate these components, using the same cutoff frequencies as before for lowfrequency (<1 kHz) and highfrequency (>2.5 kHz) EFP.
The same simplified model used in Figure 3 was fit to the data (example shown in Figure 4) by performing a nonlinear least squares optimization. The model considered only the average membrane potential across the fibers, and we calculated the membrane currents based on the density of fibers instead of simulating individual fibers. The model also discarded the radial extent of the bundle, treating it as a line; see Materials and methods for more details on the model. Free parameters to be fit were (1) the number of fibers at the depth of each recording location, (2) the average spatial derivative of the membrane potential over time in the fibers at the location next to the topmost, dorsal electrode, (3) the axonal conduction velocity, and (4) the average distance between the axon bundle and electrode array.
The resulting EFP responses and the model fit are depicted in Figure 5. Figure 5A shows the inferred average over trials of the deviation of the membrane potential from the resting potential in response to the stimulus, at a location in the axon next to the first electrode (penetration depth 1550 µm), obtained from the fit. The inferred voltage is composed of high and lowfrequency components similar to those observed in the EFP. The inferred number of fibers as a function of dorsoventral depth is shown in Figure 5B. The number (scaled by an arbitrary factor) has its maximum at the center of the electrode array, and decays steadily to both sides. This profile of the number of fibers is consistent with the known anatomy of axons in NL (Carr and Konishi, 1990; Kuokkanen et al., 2010).
The lowfrequency ($<1$ kHz, Figure 5C) responses reveal the typical polarity reversal that we predicted for an axonal terminal zone (Figures 1 and 2). The dorsoventral depth in Figure 5C and D is on the vertical axis, which corresponds to the horizontal axis in Figure 5B. The orange lines indicate the actual responses in the data.
The lowfrequency responses at the dorsal and ventral edges in Figure 5C show the same shape, but with opposite polarity, as expected for a dipolar field. Note that for a pure dipole field, the amplitude of the central responses have zero amplitude. In the data shown here, central responses show a diminished maximum amplitude, which we interpret as the contribution of higherorder (mostly quadrupole) components. The model is able to capture the behaviour of this quadrupolar component as well, with a slight underestimation of the amplitude of the peak at ventral locations. The model even captures a small oscillation in the data with period of $\approx 1$ ms in the center of the recording. Here, too, the small deviations are likely due to slightly inhomogeneous conduction velocities or nonaxonal sources.
In addition to the dipolar behavior of the lowfrequency response, we also examined the highfrequency ($>2.5$ kHz) response, shown in Figure 5D. The responses have a Gaborlike shape, as expected (Wagner et al., 2009), with maximum amplitude in the center of the recording array, at around 850 µm penetration depth. The axonal conduction velocity was calculated to be 4.0 m/s, and the distance from the bundle was 162 µm. A previously published estimate of the axonal conduction velocity in this nucleus (McColgan et al., 2014) gave a confidence bound of 0.4–6 m/s. Toward the edges (<100 µm and >1400 µm), the amplitude decays. In the central region (400–1200 µm recording depth), a systematic shift in delay can be observed, while the response appears stationary in the more dorsal and ventral electrodes. The delay increases from ventral to dorsal, which is consistent with the anatomy for contralateral stimulation.
All these aspects of the data are qualitatively reproduced by the model (Figure 5C and D, black traces). The main deviation between model and data lies in a diminished amplitude of the highfrequency oscillation modelled at the most central electrode sites (Figure 5D). Because the phase shift in the central region is mainly determined by the conduction velocity, this mismatch might be due to a variable conduction velocity in the nucleus, and the constant velocity in the model. McColgan et al. (2014) showed that different conduction velocities exist in the core and periphery of the nucleus, as predicted by Carr and Konishi (1990) from variable internode distances. A diminished amplitude in the fit could reflect an inability of the model to exactly match the phase progression. Another possible explanation is that the additional amplitude could be due to nonaxonal sources such as synaptic currents or postsynaptic spikes, which do not follow the assumptions underlying our model; see the Discussion for arguments why we expect such contributions to the EFP to be small.
When comparing the inferred membrane potential response (Figure 5A) to the measured EFP response (Figure 5C and D), the most salient difference is the dissimilar sizes of the frequency components. In the EFP, the lowfrequency component has a comparable amplitude to the highfrequency component, but in the membrane potential the lowfrequency component is much larger. This is because the EFP is related to membrane currents, which are proportional to the first and second derivatives of the membrane potential, and taking the derivative is equivalent to applying a highpass filter.
We performed the fitting procedure (example in Figure 5) for 26 recordings from 3 different owls, with monaural stimulation from both ears (implying the activation of distinct axonal populations). The average correlation coefficient for all recordings was ${R}^{2}=0.56\pm 0.15$. The correlation coefficient for the example shown in Figure 5 was 0.62.
Dipole moments of idealized axon bundles
We have shown theoretically and experimentally for specific examples of axonal projection zones and inputs how dipolar EFPs emerge. We now generalize this approach to predict the resulting dipolar EFP for arbitrary axon and stimulus configurations. Based on our cabletheory model, we analytically derived the maximal dipole moment ${p}_{\text{max}}$ for a large range of scenarios. From a given dipole moment the maximum far field potential at distance $r$ can be calculated as ${\varphi}_{\text{max}}=\frac{{p}_{\text{max}}}{4\pi {\sigma}_{e}{r}^{2}}$ where ${\sigma}_{e}$ is the extracellular conductivity.
To simplify the analytical derivation as much as possible, we assumed a Gaussian waveform for the membrane potential of a single spike, with an amplitude ${\overline{V}}_{\text{spike}}$ and a width ${\sigma}_{\text{spike}}$. The resting membrane potential was irrelevant because only the first and second derivatives of the membrane potential contribute. The axon bundle population consisted of fibers with radius $a$, axial resistance ${r}_{L}$, and conduction velocity $v$. The population was assumed to be driven with a Gaussian firingrate pulse with maximum firing rate ${\overline{\lambda}}_{\text{pulse}}$ and width ${\sigma}_{\text{pulse}}$. The distribution of the number of fibers at a given depth location was also described with a Gaussian, with width ${\sigma}_{n}$ and maximum number $\overline{n}$. This is an adequate approximation if the spikes in the incoming fibers contribute little to the dipole moment before reaching the projection zone. In this scenario, we calculated the peak dipole moment of the bundle (see Materials and methods for details) to be
Equation 1 tells us that the dipole moment is proportional to ${a}^{2}$, $\overline{n}$, ${\overline{\lambda}}_{\text{pulse}}$, ${\overline{V}}_{\text{spike}}$, and $1/{r}_{L}$. The dependence on $v$ and the widths is more complicated; the response is maximal with respect to the three (spatial) widths ${\sigma}_{n}$, $v{\sigma}_{\text{pulse}}$ and $v{\sigma}_{\text{spike}}$ when they satisfy the condition ${w}_{1}^{2}={w}_{2}^{2}+{w}_{3}^{2}$ where ${w}_{1}$ is the largest of the three terms, while ${w}_{2}$ and ${w}_{3}$ are the other two terms, regardless of order. The dipole moment is thus maximal when the widths of the spike, the pulse, and the terminal zone agree. In particular, if ${\sigma}_{n}$ (the width of the terminal zone) is the widest, then the dipole moment is maximal if ${\sigma}_{n}$ is equal to the spatial width of the overall activity in the axons, which is $v\sqrt{{\sigma}_{\text{spike}}^{2}+{\sigma}_{\text{pulse}}^{2}}$. The widths add in this way because the overall activity is the convolution of two Gaussians.
Using this formula, it is then possible to calculate the expected contributions to the EFP for different scenarios. To test the approximation in the case of the barn owl, we chose the following values: axon radius $a$ = 1 µm, conduction velocity $v$ = 4 $\frac{\text{m}}{\text{s}}$ as inferred in the previous section, axial resistivity ${r}_{L}$ = 1 $\mathrm{\Omega}\text{m}$, and extracellular conductivity ${\sigma}_{e}$ = 0.33 $\frac{\text{S}}{\text{m}}$ as used in studies of the cortex (Gold et al., 2006; Holt and Koch, 1999), anatomical and physiological parameters ${\sigma}_{n}$ = 500 µm, $\overline{n}$ = 80000, ${\overline{V}}_{\text{spike}}$ = 70 mV from (Carr and Konishi, 1990), and activation patterns for click stimulation from (Köppl, 1997a; Carr et al., 2016): ${\overline{\lambda}}_{\text{pulse}}$ = 1000 spikes/s, ${\sigma}_{\text{spike}}$ = 250 µs, ${\sigma}_{\text{pulse}}$ = 0.5 ms. This leads to a value for the dipole moment of ${p}_{\text{max}}\approx 2.5\text{\xb5A}\cdot \text{mm}$. At a distance of 750 µm, roughly the furthest distance recorded with the multielectrode array in Figure 4 and Figure 5, this dipole moment corresponded to a field potential of 1.1 mV, consistent with the order of magnitude of the responses in our experiments (Figures 4 and 5).
Dipole sources are also to be expected to make up the majority of the electrical signals recorded at the scalp (Nunez and Srinivasan, 2006). One such signal is the auditory brainstem response (ABR), which is recorded at the scalp in response to auditory stimulation with clicks or chirps (Riedel and Kollmeier, 2002). An amplitude of about 10 µV of the ABR in the barn owl has recently been reported by PalancaCastan et al. (2016). We calculated the contribution expected from an axon bundle with the same characteristics as described before at 2 cm from NL, aiming to estimate the contribution to the ABR. Multiplying by a factor of 2 to account for the fact that there is an NL in each hemisphere, the predicted contribution was 3.1 µV, which is of the same order of magnitude as the value reported in the experiments.
To estimate the lowfrequency dipole moment of NL from our multielectrode recordings, it is sufficient to use CSD analysis in one dimension, i.e. $\frac{{\partial}^{2}}{\partial {z}^{2}}\varphi (z)=\frac{1}{{\sigma}_{e}}i(z)$ and a simple sinusoidal approximation of the voltage within NL: $\varphi (z)={\varphi}_{0}\mathrm{sin}(2\pi z/L)$ for $L/2<z<L/2$ and $\varphi (z)=0$ otherwise, where ${\varphi}_{0}\approx 0.5$ mV is the amplitude, $L\approx 2$ mm is the spatial wavelength, and $z$ is the depth in NL with $z=0$ being in the center. To convert the current density $i$ into a current, we approximate the NL volume that contributes to the dipole as ${V}_{\text{NL}}\approx 6$ mm${}^{3}$ (Kuokkanen et al., 2010). We assume that the current is homogeneously distributed in the directions perpendicular to $z$. Using the definition of a dipole, ${p}_{max}:=\int \mathrm{d}Vi(z)z$, we can integrate over the dimensions perpendicular to $z$ and obtain ${p}_{max}=\frac{{V}_{\text{NL}}}{L}{\int}_{L/2}^{L/2}\mathrm{d}zi(z)z$. Substituting $i(z)$ and solving the integral, we find the maximum dipole moment to be ${p}_{max}=2\pi {V}_{\text{NL}}{\sigma}_{e}{\varphi}_{0}/L\approx 3\text{\xb5A}\cdot \text{mm}$, which is consistent with our previous estimates.
As a second example, we considered thalamocortical projections, for which (Swadlow and Gusev, 2000) reported amplitudes of extracellular spikerelated potentials, called axon terminal potentials, at various locations; for example, at 400 µm from the center of the dipole, they reported an amplitude of the response of $\approx $ 1 µV. Individual thalamocortical axons are thin and have large and highly branched projection zones (Feldmeyer, 2012), so we estimated ${\sigma}_{n}$ = 250 µm, $\overline{n}$ = 30, and $a$ = 1 µm. We assumed a jitter ${\sigma}_{\text{pulse}}$ = 125 µs in the arrival time instead of a true activity pulse, and we normalized the pulse to have area 1 because we were considering a spike triggered average. The conduction velocity has been reported as $v$ = 8.5 m/s (Simons et al., 2007). Leaving all other values as in the previous approximation, we arrived at a dipole moment of $p}_{\text{max}}\approx 1.5\phantom{\rule{thinmathspace}{0ex}}\text{\xb5A}\cdot \text{\xb5m$, yielding an extracellular spike amplitude of $\approx $ 2.3 µV at the distance of 400 µm, which is of the same order of magnitude as the value ($\approx $ 1 µV) reported by Swadlow and Gusev (2000).
In cases in which the jitter ${\sigma}_{\text{pulse}}$ is longer, the dipole moment is lower. For example, for pulses evoked by a visual stimulus, the pulse durations can exceed 10 ms (Mitzdorf, 1985; Schroeder et al., 1991; Schroeder et al., 1998; Self et al., 2013). Using the same parameters as for the thalamocortical projection employed before, but increasing the number of fibers by a factor of 100, increasing the width of the pulse to 10 ms, and increasing the maximal firing rate to 10 Hz, we found that the value of the dipole moment was ${p}_{\text{max}}\approx 0.018\text{\xb5A}\cdot \text{\xb5m}$, which is two orders of magnitude smaller than in the case of the brief pulse discussed above. However, when we further reduced the conduction velocity to 0.4 m/s, the same 10 ms pulse produced a dipole moment of $0.39\text{\xb5A}\cdot \text{\xb5m}$. Such low conduction velocities can, for example, be found in corticocortical projections (Swadlow, 1989).
To summarize, Equation 1 quantitatively predicts the contribution of axonal projection zones to the far field EFP, and this prediction matched experimental values in several cases.
Discussion
Numerical simulations, analytical calculations, and experimental data allow us to show how axonal fiber bundles may contribute to the EFP, and explain how the contributions are shaped by axonal morphology. There are three principal effects of axon bundle structure on the EFP. First, the lowfrequency components of the EFP are governed by the densities of bifurcations and terminations and can have a dipolar structure (Figure 1 and Figure 2A–C). Second, the highfrequency components are governed by the local number of fibers (Figure 2D–F). Third, membrane potentials that change on the same spatial scale as an axonal projection zone through which they propagate generate strong dipole moments in the EFP response (Figure 3). At the temporal frequencies that correspond to wavelengths of the size of the projection zone, this leads to dipolar EFP components that are not negligible and exceed the reach of the presumed quadrupolar nature of axonal EFPs.
Relevance to the interpretation of electrophysiological recordings
Our findings relate to the interpretation of a wide range of electrophysiological data in general, and to the estimation of current sources in particular. When performing a typical current source density (CSD) analysis, the local number of fibers cannot be disentangled from membrane current density (Nicholson, 1973; Potworowski et al., 2012). In CSD analysis, the membrane current densities can vary independently with time and location. In contrast, in the case of axon fiber bundles, the situation is different: the number of fibers is variable in space, in particular in the terminal zone, but the current sources at different locations are highly correlated because they are caused by propagating action potentials. In the case presented here (Figure 5) where axonal action potentials dominate the EFP, it was possible to recover the (normalized) fiber densities and average membrane potentials from the recordings.
Beyond recovering actual fiber densities and membrane potentials, our approach enables the interpretation of CSD results in the presence of axon fiber bundles. For example, the sink and source distribution found in classical CSD analysis of axon bundles (Mitzdorf and Singer, 1978; Mitzdorf, 1985; Mitzdorf and Singer, 1977) shows a dipolar structure in terminal zones, but a conclusive explanation of their origin was not given. (Tenke et al., 1993) studied the dipole at an axonal terminal zone in the macaque striate cortex for a fixed point in time, attributing the sinks to the depolarized axon endings, and the sources to the return currents distributed along the axons, while not taking account of additional currents flowing at bifurcations. Our modeling approach provides a novel way of interpreting these findings in terms of actively propagated action potentials in a fiber bundle.
As an example case for a fiber bundle, we recorded from the barn owl nucleus laminaris. Figure 4 and Figure 5 showed that the low (<1 kHz) and highfrequency (>2.5 kHz) components exhibit qualitatively different behaviours as a function of recording location relative to the terminal zone. The lowfrequency component is a largely stationary phenomenon, while the fine structure of the highfrequency component shifts gradually in space as a function of the axonal conduction velocity (Figure 5, see also (Carr et al., 2015)). Lowfrequency components have a strong dipole moment, meaning that they contribute to the farfield EFP. Due to the difference in reach, the highfrequency component is most suitable for the study of local phenomena while the lowfrequency component bears information about locations more distant from the recording site (Figure 3), consistent with findings on nonaxonal EFPs (Pettersen and Einevoll, 2008; Łęski et al., 2013).
Note that the lowpass (<1 kHz) filtered EFP is calculated in a similar way to the LFP, with the exception that the cutoff frequencies used to separate the low and highpass filtered EFP are relatively high compared to those used in cortical or hippocampal studies to separate LFP and MUA. We applied these high cutoffs because our modeling and experiments were performed in the auditory brainstem of the barn owl, which operates on very short time scales and, consequently, higher frequencies. We expect other systems operating on slower time scales to have lower optimal cutoff frequencies separating the components. Equation 1 indicates how the different spatial and temporal system properties relate to each other to generate a dipole moment.
Dipolar fields are essential for the generation of electrical field potentials at greater distances from the brain. The most prominent of these is the EEG, which is commonly attributed to the dipolar contributions of pyramidal cells (Nunez and Srinivasan, 2006). As originally suggested by Tenke et al. (1993), we propose that axonal contributions might also be relevant in the analysis of these fields. This is particularly true for the auditory brainstem response (ABR), which is closely related to the EEG and involves brain structures that display high degrees of synchrony as well as axonal organization, and are thus ideal candidates for the generation of axonal field potentials visible at long ranges. This would in turn have implications for the interpretation of the ABR in clinical contexts.
The ABR amplitude of the barn owl has been reported to be on the order of 10 µV (PalancaCastan et al., 2016) while we estimated a contribution of about 3 µV amplitude from the incoming axons in NL alone. This estimate of 30% axonal contribution to the ABR suggests that there may indeed be measurable components due to axons in the ABR, in particular, and the EEG, in general. However, this estimate is crude because it did not take into account the anatomy of the skull except for its size. Future studies based on a more detailed skull model and paired recordings of ABR and EFP should improve our understanding of axonal contributions to the ABR.
We have shown that the EFP in the barn owl NL is consistent with a model of axonal sources. We believe synaptic contributions to be small in this case for the following reasons: The somatic membrane potentials due to synaptic currents are much smaller than the impact of postsynaptic spikes (Ashida et al., 2007; Funabiki et al., 2011). Since postsynaptic spikes contribute little to the EFP (Kuokkanen et al., 2010), we suspect that the synaptic contributions to the EFP are also small. Furthermore, synaptic EFP contributions would require a spatial separation of currents, which is not possible to achieve in NL because of the symmetrical arrangement of synapses on the spherical NL cell bodies (Carr and Konishi, 1990), meaning that synaptic sources can also not explain a dipolar EFP, and are thus expected to contribute little to the EFP.
Dipolar EFPs in other animals and brain regions
It is interesting to note that a similar dipolelike reversal of polarity as shown here for the barn owl NL has been reported in the chicken NL (Schwarz, 1992) as well as in the mammalian analog to NL, the medial superior olive (MSO) (Mc Laughlin et al., 2010). The phase reversal in this case was modeled based on the assumption that the postsynaptic NL and MSO dendrites with their bipolar morphology generate the dipolar EFP response (Mc Laughlin et al., 2010; Goldwyn et al., 2014). However, in the owl this dipolar morphology of neurons is largely absent (Carr and Konishi, 1990), making dendritic sources unlikely. This differential morphology suggests that similar dipolar field potentials in owls and mammals emerge from different physiological substrates. Such a convergence might point towards an evolutionary pressure favoring a bipolar EFP structure in coincidence detection systems, and indeed, (Goldwyn and Rinzel, 2016) have proposed a model in which this extracellular potential enhances the function of coincidence detectors through a form of ephaptic interaction. Their approach centers on dendrites and is not directly transferable, but it seems possible that a similar mechanism might arise in the barn owl NL based on axonal EFPs.
The key assumption underlying our modeling of axonal geometries is that there exists a preferential direction of the axon arbor. In many structures this is the case, for example in the parts of the auditory brainstem we studied here. In other brain regions, this tendency is not as prominent, with a spectrum existing between completely directed and undirected growth. More undirected growth would lead to a more diffuse response in the EFP, and eventually to a cancellation of the dipolar field potentials. Cuntz et al. (2010) and Budd et al. (2010) studied the principles underlying the growth patterns of axons and found that the degree of direction in the growth of an axon depends on the balance struck between conduction delay and wiring cost. Optimizing for minimal conduction time leads to highly directed structures while optimizing wiring cost leads to more tortuous, undirected growth. This insight suggests that directed structures  and thus also strong, dipolar EFPs due to axons  may be more prevalent in systems which require high temporal precision in the information processing. This requirement for high temporal precision also aligns well with our model prediction: the dipole moment (Equation 1) is maximal when the spatial spread of the activation is equal to the size of the projection zone, favoring short activation times (<1 ms) for typical projection zone sizes (1 mm) and conduction velocities (1 m/s).
Relationships to more detailed biophysical models
In the systems we were aiming to describe with our model, for example NL and thalamocortical projections, synaptic boutons are typically small, and we did not model them explicitly. In other systems such as the neuromuscular junction, the synaptic structure can be very large when compared to the axon bundle (Harris and Ribchester, 1979; Katz, 1961; Katz and Miledi, 1965). Such a large junction with an overall length of up to 1 mm was modeled by Gydikov and Trayanova (1986). They found a significant effect of this structure on the EFP. The single flaring and tapering neuromuscular junction in their model had a comparable effect as the entire projection zone in our model, with the flaring causing a similar effect as the bifurcations, and the tapering taking the role of the terminations in our model. Given that synaptic boutons are several orders of magnitude smaller in NL and cortex, we do not expect a strong effect in these systems.
Membrane currents flowing in boutons were studied by Geiger and Jonas (2000), who recorded from the terminals of hippocampal mossy fibers and examined calcium and potassium conductances. The potassium conductances broadened the incoming spikes in an activitydependent manner. This spike broadening is hypothesised to be mediated by slow inactivation of the potassium channels and takes place on a timescale of >100 ms, and is thus not relevant to the present study. Spike broadening could be captured in our model by incorporating in Equation 1 a ${\sigma}_{\text{spike}}$ that is variable in time.
The calcium currents reported by Geiger and Jonas (2000) were further quantified by Alle et al. (2009). Calcium currents were temporally overlapping and much smaller in amplitude than sodium and potassium currents. We therefore neglected calcium currents in our model.
Modeling the myelinated compartments, we assumed that they are purely passive and strongly insulated from the extracellular space. However, myelinated compartments do in fact express active conductances, in particular in the paranodal and juxtaparanodal region (Chiu and Ritchie, 1981; Waxman and Ritchie, 1985). Including such a detailed distribution of ion channels in our model could lead to a different shape of the waveform and the spectrum of the EFP of an action potential, possibly similar to the effect described by Ness et al. (2016) for active conductances on dendrites. The conclusions drawn by our model are, however, independent of the precise active conductances and the distribution of myelinated and active segments along the axons because our results rely only on the gross waveform of propagating action potentials, but not on finer details. Active conductances and capacitive currents in the myelinated segments could affect the shape of the action potential waveform, but do not affect our conclusion about the spatial scaling behaviour of the EFP.
Because of the weak dependence of our results on the gross extracellular spike waveform, our analytical model does not include any intrinsic lowpass filtering as can be derived, for example, for dendritic models (Lindén et al., 2010; Einevoll et al., 2013 ; for reviews see Buzsáki et al., 2012). The effective additional currents flowing at bifurcations and terminations are, however, lowfrequency contributions to the overall membrane currents in our model. Extending our model to treat these currents separately might show whether axons could contribute to the observed $1/f$ scaling of the spectrum of the EFP (Pritchard, 1992).
Conclusion
Axonal projections can contribute substantially to EFPs. Our results quantitatively show how the anatomy of axon terminal zones and the activity in axons determine its frequencyspecific farfield contribution to the EFP.
Materials and methods
Experimental recordings
Request a detailed protocolThe experiments were conducted at the Department of Biology of the University of Maryland. Data was collected from three barn owls (Tyto furcata pratincola). Procedures conformed to NIH Guidelines for Animal Research and were approved by the Animal Care and Use Committee of the University of Maryland. Anaesthesia was induced prior to each experiment by intramuscular injection of a total of $810$ ml/kg of $20\%$ urethane divided into three to four injections over the course of $3$ hours. Body temperature was maintained at ${39}^{\circ}$C by a feedbackcontrolled heating blanket.
Recordings were made in a soundattenuating chamber (Acoustic Systems Inc., Austin, TX, USA). Tungsten electrodes with impedances between $2$ and $20$ M$\mathrm{\Omega}$ were used to find suitable recording locations in nucleus laminaris (NL). Once NL had been located, the tungsten electrode was retracted and replaced with a 32 channel multielectrode array (A1 × 32–15 mm50413A32, Neuronexus, Ann Arbor, MI, USA). The multielectrode array was lowered using a microdrive (MP225, Sutter Instruments Co., Novato, CA, USA) during continuous presentation of a whitenoise burst stimulus until visual inspection of the waveform showed that NL was at the center of the array. A grounded silver chloride pellet, placed under the animal’s skin around the incision, served as the animal ground electrode. Electrode signals were amplified by a headstage (HS36, Neuralynx Inc., Tucson, AZ, USA). An adapter (ADPTHS36N2T32A, Neuralynx Inc.) was used between the electrode and the headstage. A further adapter (ADPTHS36ERP27, Neuralynx Inc.) was used between the headstage and the control panel in order to map all 32 channels to the amplifiers. Preamplified electrode signals were passed to the control panel (ERP27, Neuralynx Inc.), then to four 8channel amplifiers (Lynx8, Neuralynx Inc.), and then to an analoguetodigital converter (Cheetah Digital Interface, Neuralynx Inc.) connected to a personal computer running Cheetah5 software (Neuralynx Inc.) where the responses were stored for offline analysis.
Acoustic stimuli were digitally generated by a custommade Matlab (MathWorks, Natick, MA, USA; RRID:SCR_001622) script (McColgan and Liu, 2017) driving a signalprocessing board (RX6, TuckerDavis Technologies, Alachua, FL, USA) at a sampling rate of 195.3125 kHz. The sound stimuli were attenuated using a programmable attenuator (PA5, TuckerDavis Technologies). Click stimuli were generated as a single halfwave of a 5 kHz sine tone. Miniature earphones were inserted into the owl’s ear canals and fixed to a headplate. Acoustic stimuli were fed to these earphones. Stimulus delivery was triggered by National Instruments equipment (NI USB6259 and BNC2090A, National Instruments Inc, Austin, TX, USA), and stimulus presentation times were recorded along with the responses. Trigger pulses were configured in MatLab through Ephus software (Vidrio Technologies LLC, Ashburn, VA, USA). Responses were averaged over 10 repetitions.
Multicompartment model of axons
Request a detailed protocolWe modeled axons using NEURON (Hines and Carnevale, 1997; Hines et al., 2009) and extended previous work by Simon et al. (1999) and Kuba and Ohmori (2009), which included the high and lowthreshold potassium channels used by Rathouz and Trussell (1998). The axon was modeled as a sequence of active nodes and passive myelinated segments. The nodes and myelinated segments had lengths of 2 µm and 75 µm, respectively. We used the model of a nucleus magnocellularis (NM) axon provided by Simon et al. (1999) in ModelDB (Hines et al. (2004), Accession number: 153998). In order to ensure robust spike propagation at the bifurcations, some of the model the parameters were modified. The values of the properties that were modified from those provided by Simon et al. (1999) are shown in Table 1. In addition, the Q10 values were set to 3, and the temperature was set to ${40}^{\circ}\text{C}$ as done by Kuba and Ohmori (2009). The ratio of leak conductance and capacitance between node and myelin was changed from 80 as used by Simon et al. (1999) to 1000 (Koch, 2004). We removed the HodgkinHuxley type potassium conductivities included by Simon et al. (1999) (which are based on data from the squid axon) from the simulations, and increased the KLVA and KHVA conductances (which are based on physiological data from the auditory brainstem) to compensate. In order to avoid nodes lining up in axial direction, the very first myelinated segment had a length drawn from a uniform distribution between 0 and 75 µm.
We included branching axons in our simulations. Branches were generated by connecting the incoming passive segment to one end of a node, and the two outgoing passive segments to the other end of the node, and then continuing the alternation of active and passive segments in each resulting branch. In Figure 1A–D, where the positions of bifurcations or terminations of axons were fixed, the last node was placed in the axon as usual, and the last myelinated segment was shortened in order to obtain the total length before the bifurcation or termination.
To approximate the infinitely long axons, we added segments before and after the shown portions of the axon. We chose the total length of the additional segments by incrementally increasing the length segmentbysegment until there was no visual difference between each successive lengthening of the axon, arriving at a length of 3 mm at each end.
To evoke an action potential at a designated time, a special conductance was temporarily activated in the first node of Ranvier. The conductance had a reversal potential of 0 mV, a maximal amplitude of 0.05 µS, and a time course described by an alphafunction with time constant 0.01 ms. Soma and axon initial segment were not modeled explicitly. This conductance resembled a synaptic conductance, except for the very short time constant and the lack of a somatic or dendritic compartment on which synapses typically impinge.
For the simplified axon geometries used in Figure 1, the branching pattern was fixed as described in the caption, with the exception of the axial positions of the branching points in Figure 1E: a random offset between branching points was drawn from a gamma distribution with mean 400 µm and standard deviation 300 µm. The initial branching point for each axon was offset from the original location by a distance drawn from a Gaussian distribution with mean zero and a width of 300 µm. This was done to smooth out the effects of individual branchings or terminations.
For the axons in Figure 2, branching patterns were generated procedurally, starting with a root segment. In order to avoid artifacts from the stimulus and to simulate a long fiber tract prior to the terminal zone, a sequence of 10 active and passive segments without bifurcations was assumed before the terminal zone (770 µm total length). To this root, segments were appended iteratively. Before adding a segment, a decision whether to branch or terminate an axon was drawn from a probability distribution that was dependent on the axial position of the end of the previous segment. These probability distributions were modeled as logistic functions with the parameters adjusted to roughly match the numbers of branchings and terminations shown by Carr and Konishi (1990). Thus, an initial phase dominated by bifurcations was followed by a phase dominated by terminations, with the probability of termination reaching 100% at the end of the terminal zone. The distribution of bifurcations had its maximum at axial location $z=0$ with a standard deviation of 200 µm. The distribution of terminations had its maximum at $z=500$ µm, with a standard deviation of 100 µm. The branching angle had an average of ${20}^{\circ}$, with a standard deviation of ${5}^{\circ}$. At branching points, the plane containing the branches had a uniform random orientation, resulting in a 3dimensional structure of the axon bundle.
In all cases except for the simulations shown in Figure 2, the radial position of all node and myelin compartments was set to zero, meaning they were placed on a straight line extending axially.
Numerical simulations of action potentials propagating along axons yielded the membrane currents from which we calculated extracellular fields. This procedure is described in detail by Gold et al. (2006). Briefly, the extracellular medium is assumed to be a homogeneous volume conductor with conductivity ${\sigma}_{e}$, and a quasistatic approximation of the electrical field potential $\varphi $ is made. The extracellular potential $\varphi (\mathbf{\mathbf{r}},t)$ at location $\mathbf{\mathbf{r}}$ and time $t$ due to a membrane current density distribution $i(\mathbf{\mathbf{r}},t)$ is then governed by the equation $\mathrm{\Delta}\varphi (\mathbf{\mathbf{r}},t)=\frac{1}{{\sigma}_{e}}i(\mathbf{\mathbf{r}},t)$, with $\mathrm{\Delta}$ denoting the Laplace operator Nunez and Srinivasan (2006). If the currents $i$ are constrained to a volume $W$, this equation has the solution
Since the majority of the current flows through the small nodes of Ranvier in myelinated axons, we used the pointsource approximation for all compartments, and subdivided the myelinated segments into 10 isopotential sections each; we did not use the linesource approximation (Holt and Koch, 1999).
Analysis of the resulting extracellular field potential (EFP) included filtering. All filtering was performed with thirdorder Butterworth filters. The multi unit activity (MUA) was calculated by highpass filtering the signal with a cutoff of 2500 Hz, setting all samples with negative values to zero, and then lowpass filtering the resulting response with a cutoff of 500 Hz. The lowpass filtered EFP was calculated with a cutoff of 1000 Hz. To exclude influences of spectral leakage on our results, we also performed the same analysis with 20thorder Butterworth filters and found qualitatively identical results. The specific filtering (3rd or 20th order Butterworth) did not affect our conclusions.
The code for these simulations is available at https://github.com/phreeza/pyLaminaris (McColgan and Liu, 2017). A copy is archived at https://github.com/elifesciencespublications/pyLaminaris.
Meanfield model of an axon bundle
Request a detailed protocolTo better understand the processes leading to the complex spatiotemporal patterns of extracellular fields, we devised an analytically tractable model of axon bundles. We defined the spatial dimension in cylindrical coordinates $\mathbf{\mathbf{r}}=(\rho ,\phi ,z)$, where $\rho $ was the radial distance from the cylindrical axis, $\phi $ the angle of azimuth, and $z$ the axial distance along the cylinder axis. We considered a simple model axon bundle that extended infinitely in the axial $z$direction at $\rho =0$. The bundle had a variable number of fibers along the $z$ axis, denoted by $n(z)$, each of which cylindrical with an identical radius $a$. This meant that the total crosssectional area $A$ of the bundle at a given depth $z$ was given by $A(z)=\pi {a}^{2}n(z)$. We assumed the axons to be perfect transmission lines, meaning that the action potential is a traveling wave with velocity $v$ along the axon. In particular, we neglected delays and distortions that can be induced when an axon branches or terminates. In this case, we could assume that the membrane voltage was the same in each fiber for a given $z$ coordinate. From linear cable theory (e.g. Jack et al., 1975), we then obtained the following expression for the total transmembrane current per unit length $I(z,t)$ from a given membrane potential $V(z,t)$:
We next calculated the corresponding extracellular field potential $\varphi (\mathbf{\mathbf{r}},t)$ of a given membrane potential waveform $V(z,t)$ propagating through the axon bundle. Due to the rotational symmetry of the system and the fact that current flows only at $\rho =0$, the membrane current can be described as $i(\mathbf{\mathbf{r}},t)=I(z,t)\frac{\delta (\rho )}{\rho}$, where $\frac{\delta (\rho )}{\rho}$ is the Dirac delta distribution for a line at $\rho =0$. Applying this membrane current to Equation 2, we obtained
which is independent of the angle $\phi $.
To derive the frequency responses in Figure 3, we simulated the membrane potentials as sinusoids, that is, $V(z,t)=\mathrm{sin}\left(2\pi f\cdot \left(ztv\right)\right)$, with varying frequencies $f$ between 100 Hz and 5 kHz and calculated the standard deviation of the response for each frequency individually. The amplitude of the membrane potential $V$ was the same for all frequencies.
To derive the dipole moment of a simplified projection zone, we considered an axon bundle in which identical spikes with the waveform ${V}_{\text{spike}}(z,t)$ propagate as traveling waves with a velocity $v$ in positive $z$ direction: ${V}_{\text{spike}}(z,t)={V}_{\text{spike}}(ztv,0)$. If each of the fibers is stimulated with an inhomogeneous Poisson process, with the driving rate $\lambda (t)$ shared among all axons, the average membrane potential across fibers will be $V(z,t)={V}_{\text{spike}}(z,t)\ast \lambda (t)$, where $\ast $ denotes the convolution with respect to the time $t$. Plugging this into Equation 5, we obtained
Analytical solution of the meanfield model of an axon bundle
Request a detailed protocolAssuming Gaussian shapes for the firingrate pulse $\lambda (t)={\overline{\lambda}}_{\text{pulse}}\mathrm{exp}\left(\frac{{t}^{2}}{2{\sigma}_{\text{pulse}}^{2}}\right)$, the spike ${V}_{\text{spike}}(z,t)={V}_{\text{spike}}(ztv,0)={\overline{V}}_{\text{spike}}\mathrm{exp}\left(\frac{{(ztv)}^{2}}{2{\sigma}_{\text{spike}}^{2}}\right)$, and the fiber bundle projection zone $n(z)=\overline{n}\mathrm{exp}\left(\frac{{z}^{2}}{2{\sigma}_{\text{n}}^{2}}\right)$, we were able to take advantage of the fact that the product and the convolution of two Gaussians are again Gaussian, and obtained
The dipole moment $p(t)$ is defined as
into which we can enter Equation 9 to obtain
The dipole moment has its peak amplitude at ${t}_{\text{max}}=\sqrt{{\sigma}_{n}^{2}+{v}^{2}\left({\sigma}_{\text{pulse}}^{2}+{\sigma}_{\text{spike}}^{2}\right)}/v$, and takes the value
which is presented as Equation 1 in the Results section.
Model fitting to experimental data
Request a detailed protocolIn order to relate the model to experimentally obtained data as shown in Figure 5, we performed a nonlinear least squares fit, minimizing the mean squared error $\u03f5$ between the measured potential ${\varphi}_{\text{measured}}$ and the model prediction ${\varphi}_{\text{model}}$ in Equations 6 and 5 for $N=32$ measurement locations ${z}_{n}$ ($n=1,\mathrm{\dots},N$) and $M=600$ time points ${t}_{m}$ ($m=1,\mathrm{\dots},M$): $\u03f5=\frac{1}{NM}\sum _{n=1}^{N}\sum _{m=1}^{M}{\left[{\varphi}_{\text{measured}}({z}_{n},{t}_{m}){\varphi}_{\text{model}}({z}_{n},{t}_{m})\right]}^{2}$. The separation between electrodes was given by the electrode layout as 50 µm. The time between sampling points was 5.12 µs. We achieved the minimization of the error $\u03f5$ using the ‘optimize.minimize’ routine provided by the SciPy package (Jones et al., 2001). The free parameters to be determined by the optimization routine were the distance $\rho $, the velocity $v$, the number of fibers per unit length $n({z}_{n})$ for each measurement location ${z}_{n}$, and the spatial derivative of the average membrane potential $\frac{\text{d}}{\text{d}z}V({z}_{1},{t}_{m})$ at electrode location ${z}_{1}$ for each time point ${t}_{m}$. We fit the first derivative of the membrane potential in order to better capture the lowfrequency components that we found in Figure 1E, and because the membrane potential appears only as the derivative in the model. The derivative of the membrane voltage at the other locations than ${z}_{1}$ was then determined by the travelingwave assumption: $\frac{\text{d}}{\text{d}z}V({z}_{n},{t}_{m})=\frac{\text{d}}{\text{d}z}V({z}_{1},{t}_{m}\frac{{z}_{n}{z}_{1}}{v})$, using a linear interpolation between timepoints. The model assumption of a single line of axons with electrodes at a fixed distance is a simplification of a threedimensional axon tree where the fibers are distributed at various distances. The distance parameter $\rho $ in Equation 6 can be interpreted as an average distance in this simplification.
To aid the convergence of the fit algorithm, an initial guess for the number of fibers $n({z}_{n})$ was set by hand. Initializing the guess to a constant or a fully random number of fibers resulted in a failure to converge. However, different Gaussianlike initial guesses converged to a single solution, meaning that the specific initial guess did not alter the final fit result. Initializing the membrane voltage with different normally distributed values did not affect the outcome of the fit. The results shown in Figure 5 were obtained with an initial guess of a Gaussian with amplitude 12, centered at penetration depth 725 µm with standard deviation 400 µm.
Because of the linearity of Equations 26 both in the current $I$ and the membrane potential $V$, inferring the membrane voltage $V$ from the average over trials of the extracellular potential $\varphi $ produces the average membrane voltage $V$. This in turn is the membrane voltage response of a single spike convolved with the peristimulus time histogram (PSTH).
References

Simple diagrammatic rules for solving dendritic cable problemsPhysica A: Statistical Mechanics and its Applications 185:343–356.https://doi.org/10.1016/03784371(92)904745

Cell type and activitydependent extracellular correlates of intracellular spikingJournal of Neurophysiology 114:608–623.https://doi.org/10.1152/jn.00628.2014

Passive soma facilitates submillisecond coincidence detection in the owl's auditory systemJournal of Neurophysiology 97:2267–2282.https://doi.org/10.1152/jn.00399.2006

Crossfrequency phasephase coupling between θ and γ oscillations in the hippocampusJournal of Neuroscience 32:423–435.https://doi.org/10.1523/JNEUROSCI.412211.2012

Ultrarapid axonaxon ephaptic inhibition of cerebellar Purkinje cells by the pinceauNature Neuroscience 17:289–295.https://doi.org/10.1038/nn.3624

Neocortical axon arbors tradeoff material and conduction delay conservationPLoS Computational Biology 6:e1000711.https://doi.org/10.1371/journal.pcbi.1000711

The origin of extracellular fields and currentsEEG, ECoG, LFP and spikesNature Reviews Neuroscience 13:407–420.https://doi.org/10.1038/nrn3241

BookThe role of conduction delay in creating sensitivity to interaural time differencesIn: van Dijk P, Başkent D, Gaudrain E, de Kleine E, Wagner A, Lanting C, editors. Physiology, Psychoacoustics and Cognition in Normal and Impaired Hearing, Vol. 894. Advances in Experimental Medicine and Biology Springer International Publishing. pp. 189–196.

A circuit for detection of interaural time differences in the brain stem of the barn owlJournal of Neuroscience 10:3227–3246.

Maps of interaural delay in the owl's nucleus laminarisJournal of Neurophysiology 114:1862–1873.https://doi.org/10.1152/jn.00644.2015

One rule to grow them all: a general theory of neuronal branching and its practical applicationPLoS Computational Biology 6:e1000877.https://doi.org/10.1371/journal.pcbi.1000877

Modelling and analysis of local field potentials for studying the function of cortical circuitsNature Reviews Neuroscience 14:770–785.https://doi.org/10.1038/nrn3599

Excitatory neuronal connectivity in the barrel cortexFrontiers in Neuroanatomy 6:24.https://doi.org/10.3389/fnana.2012.00024

Cytoarchitectonic and dynamic origins of giant positive local field potentials in the dentate gyrusJournal of Neuroscience 33:15518–15532.https://doi.org/10.1523/JNEUROSCI.033813.2013

Computation of interaural time difference in the owl's coincidence detector neuronsJournal of Neuroscience 31:15245–15256.https://doi.org/10.1523/JNEUROSCI.212711.2011

On the origin of the extracellular action potential waveform: A modeling studyJournal of Neurophysiology 95:3113–3128.https://doi.org/10.1152/jn.00979.2005

A model of the medial superior olive explains spatiotemporal features of local field potentialsJournal of Neuroscience 34:11705–11722.https://doi.org/10.1523/JNEUROSCI.017514.2014

Cell recognition during neuronal developmentScience 225:1271–1279.https://doi.org/10.1126/science.6474176

Influence of the muscle fibre end geometry on the extracellular potentialsBiological Cybernetics 54:1–8.https://doi.org/10.1007/BF00337110

Extracellular potentials of single active muscle fibres: effects of finite fibre lengthBiological Cybernetics 53:363–372.https://doi.org/10.1007/BF00318202

Models of axon guidance and bundling during developmentProceedings of the Royal Society B: Biological Sciences 266:2231–2238.https://doi.org/10.1098/rspb.1999.0913

The NEURON simulation environmentNeural Computation 9:1179–1209.https://doi.org/10.1162/neco.1997.9.6.1179

ModelDB: A Database to Support Computational NeuroscienceJournal of Computational Neuroscience 17:7–11.https://doi.org/10.1023/B:JCNS.0000023869.22017.2e

Electrical interactions via the extracellular potential near cell bodiesJournal of Computational Neuroscience 6:169–184.https://doi.org/10.1023/A:1008832702585

Propagation of electric activity in motor nerve terminalsProceedings of the Royal Society B: Biological Sciences 161:453–482.https://doi.org/10.1098/rspb.1965.0015

The Terminations of the Afferent Nerve Fibre in the Muscle Spindle of the FrogPhilosophical Transactions of the Royal Society B: Biological Sciences 243:221–240.https://doi.org/10.1098/rstb.1961.0001

BookBiophysics of Computation: Information Processing in Single NeuronsUnited States: Oxford University Press.

The owl's cochlear nuclei process different sound localization cuesThe Journal of the Acoustical Society of America 78:360–364.https://doi.org/10.1121/1.392499

Frequency tuning and spontaneous activity in the auditory nerve and cochlear nucleus magnocellularis of the barn owl Tyto albaJournal of neurophysiology 77:364–377.

Phase locking to high frequencies in the auditory nerve and cochlear nucleus magnocellularis of the barn owl, Tyto albaJournal of Neuroscience 17:3312–3321.

Linear summation in the barn owl's brainstem underlies responses to interaural time differencesJournal of Neurophysiology 110:117–130.https://doi.org/10.1152/jn.00410.2012

On the origin of the extracellular field potential in the nucleus laminaris of the barn owl (Tyto alba)Journal of Neurophysiology 104:2274–2290.https://doi.org/10.1152/jn.00395.2010

Frequency dependence of signal power and spatial reach of the local field potentialPLoS Computational Biology 9:e1003137.https://doi.org/10.1371/journal.pcbi.1003137

Intrinsic dendritic filtering gives lowpass power spectra of local field potentialsJournal of Computational Neuroscience 29:423–444.https://doi.org/10.1007/s1082701002454

Oscillatory dipoles as a source of phase shifts in field potentials in the mammalian auditory brainstemJournal of Neuroscience 30:13472–13487.https://doi.org/10.1523/JNEUROSCI.029410.2010

A functional circuit model of interaural time difference processingJournal of Neurophysiology 112:2850–2864.https://doi.org/10.1152/jn.00484.2014

Laminar segregation of afferents to lateral geniculate nucleus of the cat: an analysis of current source densityJournal of neurophysiology 40:1227–1244.

Current sourcedensity method and application in cat cerebral cortex: investigation of evoked potentials and EEG phenomenaPhysiological reviews 65:37–100.

Active subthreshold dendritic conductances shape the local field potentialThe Journal of Physiology 594:3809–3825.https://doi.org/10.1113/JP272022

Theoretical analysis of field potentials in anisotropic ensembles of neuronal elementsIEEE Transactions on Biomedical Engineering 20:278–288.https://doi.org/10.1109/TBME.1973.324192

The binaural interaction component in barn owl (tyto alba) presents few differences to mammalian dataJournal of the Association for Research in Otolaryngology 17:577–589.https://doi.org/10.1007/s1016201605837

Amplitude variability and extracellular lowpass filtering of neuronal spikesBiophysical Journal 94:784–802.https://doi.org/10.1529/biophysj.107.111179

Action potential sources and their volume conductor fieldsProceedings of the IEEE 65:601–611.https://doi.org/10.1109/PROC.1977.10539

Kernel current source density methodNeural Computation 24:541–575.https://doi.org/10.1162/NECO_a_00236

The brain in fractal time: 1/flike power spectrum scaling of the human electroencephalogramInternational Journal of Neuroscience 66:119–129.https://doi.org/10.3109/00207459208999796

Branching dendritic trees and motoneuron membrane resistivityExperimental Neurology 1:491–527.https://doi.org/10.1016/00144886(59)900469

Characterization of outward currents in neurons of the avian nucleus magnocellularisJournal of Neurophysiology 80:2824–2835.

Field potential signature of distinct multicellular activity patterns in the mouse hippocampusJournal of Neuroscience 30:15441–15449.https://doi.org/10.1523/JNEUROSCI.253510.2010

On highfrequency field oscillations (>100 Hz) and the spectral leakage of spiking activityJournal of Neuroscience 33:1535–1539.https://doi.org/10.1523/JNEUROSCI.421712.2013

The spiking component of oscillatory extracellular potentials in the rat hippocampusJournal of Neuroscience 32:11798–11811.https://doi.org/10.1523/JNEUROSCI.065612.2012

Can central neurons reproduce sound waveforms? An analysis of the neurophonic potential in the laminar nucleus of the chickenThe Journal of otolaryngology 21:30–38.

A dendritic model of coincidence detection in the avian brainstemNeurocomputing 2627:263–269.https://doi.org/10.1016/S09252312(99)00020X

Thalamocortical conduction times and stimulusevoked responses in the rat whiskertobarrel systemJournal of Neurophysiology 98:2842–2847.https://doi.org/10.1152/jn.00800.2007

The highly irregular firing of cortical cells is inconsistent with temporal integration of random EPSPsJournal of Neuroscience 13:334–350.

Predicting movement from multiunit activityJournal of Neuroscience 27:8387–8394.https://doi.org/10.1523/JNEUROSCI.132107.2007

Segregation of stimulus phase and intensity coding in the cochlear nucleus of the barn owlJournal of Neuroscience 4:1787–1799.

Activation of a cortical column by a thalamocortical impulseJournal of Neuroscience 22:7766–7773.

The influence of single VB thalamocortical impulses on barrel columns of rabbit somatosensory cortexJournal of Neurophysiology 83:2802–2813.

Efferent neurons and suspected interneurons in S1 vibrissa cortex of the awake rabbit: receptive fields and axonal propertiesJournal of Neurophysiology 62:288–308.

Highfrequency EEG covaries with spike burst patterns detected in cortical neuronsJournal of Neurophysiology 105:2951–2959.https://doi.org/10.1152/jn.00327.2010

Influence of spiking activity on cortical local field potentialsThe Journal of Physiology 591:5291–5303.https://doi.org/10.1113/jphysiol.2013.258228
Decision letter

Mark CW van RossumReviewing Editor; University of Edinburgh, United Kingdom
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "Dipolar extracellular potentials generated by axonal projections" for consideration by eLife. Your article has been reviewed by Timo van Kerkoerle (Reviewer #1); Gaute T Einevoll (Reviewer #2); and Rishikesh Narayanan (Reviewer #3), and the evaluation has been overseen by a Reviewing Editor and Andrew King as the Senior Editor.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Summary:
The manuscript investigates the contribution of action potentials to the local field potential, and shows how, contrary to common wisdom, axonal arbors can generate dipolar fields. The manuscript is well written and the figures are generally clear.
Essential revisions:
The revisions requested are of twofold nature: 1) more exploration of the generality of the effects, 2) more justification and exploration of the findings.
As mentioned in the manuscript an important factor in the contribution of action potentials is the synchronicity between neighboring axons; if there is too much jitter in the spikes coming in on neighboring axons the fields generated by them will cancel it. It would be of interest to quantify what the synchronicity needs to be on neighboring axons to be able to see a significant effect in the LFP. Alternatively, for a fixed level of coherence between the spikes, what would the generated field be for different frequencies? It can be judged from Figure 3 that even for bifurcations on the same axon, the frequency has to be high to still see an effect. It would therefore be relevant to see for a realistic level of coherence between the spikes what the generated axonal field would be for different frequencies. This can help to understand whether this phenomenon could play a role in (1) highfrequency ripples, as for example measured with ECoG (Worrell et al. Brain 2008), or even in (2) highfrequency gamma activity, as measured as dipoles within the cortex (van Kerkoerle et al. PNAS 2014), or with MEG on the scalp (e.g. Siegel et al. Neuron 2008).
This result will also make it easier to judge how to relate it to the previous literature (which part of the LFP can be attributed to dipoles generated by axonal arbors), and thereby how relevant the current manuscript is for a broader audience.
The auditory brainstem of the barn owl is a rather special case where the activation is very brief, in the order of a few milliseconds. The authors also refer to early papers in the visual cortex (Mitzdorf and Singer 1978; Mitzdorf, 1985), but these are special cases where electrical stimulation of the optic radiation was used, which also gives a dipolar field in the order of milliseconds as well. In contrast, the dipolar fields as evoked by a visual stimulus, are generally in the order of tens of milliseconds (e.g. Mitzdorf, 1985; Schroeder et al., 1991; Schroeder, Mehta and Givre, 1998; Self et al., 2013). It would be relevant to understand what part of a more common visual activation like that can be attributed to axonal fields.
The paper makes it conceivable that population axonal spikes contribute to dipoles which would be measurable in the LFP and potentially even to scalp EEG/MEG.
Can the authors extract the dipole and quadrupole moments from the fit?
Can the authors assess the uniqueness of the fit in their optimization procedure? Could there be several disparate solutions to this fitting procedure, possibly resulting in significant variability in the computed moments based on the specific fit that was employed?
An improved comparison is needed between the numerically comprehensive scheme and the simplified analytical model. It would be nice to illustrate how similar they are for a wellchosen example.
"All filtering was performed with thirdorder Butterworth filters […] The LFP was calculated by lowpass filtering the signal with a cutoff of 1000 Hz."
The filter used by the authors to obtain the LFP is a 3rd order Butterworth with a cutoff of 1000 Hz. The cutoff frequency is large, in comparison to most studies employing a value set at 100 Hz to 500 Hz (Buzsaki, Anastassiou and Koch, 2012; Schomberg et al., 2012; Taxidis et al., 2015). Additionally, as the order of the filter is 3, the sharpness of the filter transition from the pass band to the stop band would be shallow. This implies that frequencies beyond 1000 Hz also pass through. Additionally, in some figures, the low frequency component is arbitrarily represented at an even higher 2 kHz. EFP and LFP are alternately employed, and there is a lot of ambiguity with reference to notation and terminology. The specific choice for cutoff frequencies is of utmost importance because the authors' central claim about the contributions of (high frequency) action potentials propagating through the axon to LFPs. How dependent is the filter structure for the conclusions drawn by the authors? Would the conclusions be different if the authors employed, say a 20th order Butterworth filter with a 250 Hz cutoff?
In describing the model fit, it would be valuable to describe the most important assumptions made. Also other details of modeling need attention:
a) Location of nodes of Ranvier: In the case where there were multiple axons, were the nodes of Ranvier placed at the same spatial axis?
b) "Note that in order to avoid confounding effects, the horizontal distances between axons in Figure 1CE are for illustration only; all collaterals were simulated to lie on a straight line." Isn't this an oversimplification that artificially boosts contributions of the transmembrane currents from different axons to the EFP? What happens if this is changed to reflect an anatomical organization where horizontal distances do exist?
c) What was the structure of the axon at the bifurcations – was the surface area larger at bifurcations? What was the location of the terminal with reference to the nearest node? What would be the impact of terminating axons with a bouton with appropriate dimensions containing all ion channels (calcium, sodium and potassium especially) they are supposed to express? How would the bouton geometry (Gydikov et al., 1984 and Gydikob and Trayanova, 1986) and ion channel densities there, along with changes in each parameter mentioned here, alter the conclusions?
d) The authors assume that the myelin sheath does not contribute to the EFPs. However, given the several channels and receptors that these glial cells express, and given the active roles they play during action potential propagation, is it reasonable to assume that they don't contribute to the EFPs? What would be the contribution of the capacitative current that is consequent to the dielectric myelin between two parallel conducting media, with a highfrequency action potential propagating across? Is there any contribution of the 1/f^{n} scaling that is generally observed in LFPs, especially given the absence of intrinsic dendritic filtering in the analysis here (Nature Reviews Neuroscience reviews of Buzsaki et al., 2012, Einevoll et al., 2013)? Could the capacitive current (with reference to the myelinated compartments), however small, contribute to such filtering?
e) Please provide details of how the infinite axon was simulated.
https://doi.org/10.7554/eLife.26106.009Author response
As mentioned in the manuscript an important factor in the contribution of action potentials is the synchronicity between neighboring axons; if there is too much jitter in the spikes coming in on neighboring axons the fields generated by them will cancel it. It would be of interest to quantify what the synchronicity needs to be on neighboring axons to be able to see a significant effect in the LFP. Alternatively, for a fixed level of coherence between the spikes, what would the generated field be for different frequencies? It can be judged from Figure 3 that even for bifurcations on the same axon, the frequency has to be high to still see an effect. It would therefore be relevant to see for a realistic level of coherence between the spikes what the generated axonal field would be for different frequencies.
Synchronicity between the action potentials in neighboring axons is indeed a key factor in our manuscript. One way to quantify synchronicity is the width of a population pulse. This approach was used throughout the text and in all Figures (except Figure 3). One of our main messages of the manuscript is that the spatial width of this population pulse should match the length of the terminal zone in order to generate the maximum dipolar contribution to the EFP. The explicit dependence of the dipole moment on the synchronicity (of a pulse) is given in Equation 1, where the symbol σ_{pulse} describes the width of the pulse.
Another way to study synchronicity is axonal activity oscillating at a particular frequency, as used in Figure 3. Such an oscillating axonal activity that is travelling down an axonal bundle at a particular velocity leads to a spatial oscillation pattern along the axon. This pattern can generate a dipole moment that depends on the spatial wavelength. The dipolar component in the EFP achieves a maximum if the spatial wavelength matches the width of the terminal zone (Figures 3C and D); and the dipole moment approaches zero in the limits of very small and much larger frequencies.
For a fixed spatial wavelength, one could also modify the coherence of the oscillation, which is also called the "strength of phase locking" or the "vector strength". Even though the generated dipole naturally increases with increasing vector strength, the shapes of the curves in all graphs in Figure 3 are independent of phase locking because they were normalized or show ratios. In this way we intended to disentangle the impact of the frequency on the dipole moment from the effect of coherence.
In order "to see for a realistic level of coherence between the spikes what the generated axonal field would be for different frequencies", we would like to refer to Equation 1, which quantifies the impact of the width of pulse packets on the dipole moment.
Moreover, to better illustrate the relation of coherence and frequency in Figure 3, we added in panels A and B sketches that show the spatial wavelength of the oscillation with respect to the terminal zone. In Figures 3C and 3D we added (to the horizontal frequency axis) an axis that describes the wavelength.
This can help to understand whether this phenomenon could play a role in (1) highfrequency ripples, as for example measured with ECoG (Worrell et al. Brain 2008), or even in (2) highfrequency gamma activity, as measured as dipoles within the cortex (van Kerkoerle et al. PNAS 2014), or with MEG on the scalp (e.g. Siegel et al. Neuron 2008).
This result will also make it easier to judge how to relate it to the previous literature (which part of the LFP can be attributed to dipoles generated by axonal arbors), and thereby how relevant the current manuscript is for a broader audience.
We think that it is unlikely that locally and intrinsically generated oscillations such as hippocampal ripples or γ will contribute to local EFPs, EcOGs, or EEGs of axonal origin because the resulting spatial wavelengths are rather long. For example, for 200 Hz ripple oscillations and ~0.5 m/s conduction velocity (Schaffer collaterals) the spatial wavelength is ~2.5 mm, which is much larger than typical projection zones in the hippocampal formation. Action potentials contribute more to the EFP for higher frequencies (or shorter pulses), which are more relevant for sensory driven events, for example in the auditory pathway (measured as the auditory brainstem response, ABR).
To address this issue, we have added in the Discussion at the end of the section “Dipolar EFPs in other animals and brain regions” a sentence that summarizes the dependence of pulse width, propagation velocity, and size of the projection zone, and there is also a related new paragraph at the end of the Results. Furthermore, the improved Figure 3 could aid in understanding.
The auditory brainstem of the barn owl is a rather special case where the activation is very brief, in the order of a few milliseconds. The authors also refer to early papers in the visual cortex (Mitzdorf & Singer 1978; Mitzdorf 1985), but these are special cases where electrical stimulation of the optic radiation was used, which also gives a dipolar field in the order of milliseconds as well. In contrast, the dipolar fields as evoked by a visual stimulus, are generally in the order of tens of milliseconds (e.g. Mitzdorf, 1985; Schroeder et al., 1991; Schroeder, Mehta and Givre, 1998; Self et al., 2013). It would be relevant to understand what part of a more common visual activation like that can be attributed to axonal fields.
ABRs and auditory evoked potentials (AEPs) are brief, but widely relevant for clinical and basic research on the auditory system. We hope our findings will be directly relevant to auditory studies. For longer pulses, with durations of 10 ms or more, the expected contributions to the EFP are indeed quite small. However, the dipole moment increases with decreasing axonal conduction velocity.
To better explain what the expected contributions to the EFP of more common visual stimulations are, we added a paragraph in the Results (new paragraph before the Discussion). We are also referring back to this point in the Discussion where we now outline more explicitly the impact of axonal contributions in other systems in the section "Dipolar EFPs in other animals and brain regions".
The paper makes it conceivable that population axonal spikes contribute to dipoles which would be measurable in the LFP and potentially even to scalp EEG/MEG.
Can the authors extract the dipole and quadrupole moments from the fit?
We calculated the dipole moment from the data, as suggested, and added the value of ~3 µA mm we obtained to the manuscript in the last section of the Results.
We limited ourselves to the dipole moments. Calculating quadrupole moments directly from the currents would require an additional extensive discussion of the underlying theory. We think that such a discussion would add little information because the field potential at long distances will always be dominated by the dipolar component.
Can the authors assess the uniqueness of the fit in their optimization procedure? Could there be several disparate solutions to this fitting procedure, possibly resulting in significant variability in the computed moments based on the specific fit that was employed?
We attempted to assess the uniqueness of the fit by choosing different initial conditions to the fitting procedure. We added an explanation of this to the Materials and methods section “Model fitting to experimental data”. The combination of (1) the convergence of the fit from different starting points to similar fit results and (2) the fact that the fit results agree with known anatomy and activation statistics of the system give us the confidence that this is indeed the global optimum.
An improved comparison is needed between the numerically comprehensive scheme and the simplified analytical model. It would be nice to illustrate how similar they are for a wellchosen example.
We compared EFPs generated by (A) the numerical model as shown in Figure 2, and (B) the analytical model with a matched activation and fiber distribution. In Author response image 1 we show the heatmap of the difference ratio AB/(A+B), averaged over time. The fiber bundle lies along the horizontal axis and terminates at axial distance 16500 µm. Spikes propagate from left to right.
The relative difference at distances greater than 1 mm was < 0.3 in radial direction and < 0.05 in axial direction. The small bump that radiates from the projection zone in radial direction (at axial distance ~15500 µm) is a spurious dipole. This spurious dipole is due to a differential activation of fibers that are radially separated.
We have mentioned this consistency check in the Results in the subsection "Effects of bifurcations and terminations on distance scaling of EFPs".
The spurious dipole is also further discussed in our reply to the question about under or overestimating contributions by removing the radial dimension in the axon layout.
"All filtering was performed with thirdorder Butterworth filters […] The LFP was calculated by lowpass filtering the signal with a cutoff of 1000 Hz."
The filter used by the authors to obtain the LFP is a 3rd order Butterworth with a cutoff of 1000 Hz. The cutoff frequency is large, in comparison to most studies employing a value set at 100 Hz to 500 Hz (Buzsaki, Anastassiou and Koch, 2012; Schomberg et al., 2012; Taxidis et al., 2015). Additionally, as the order of the filter is 3, the sharpness of the filter transition from the pass band to the stop band would be shallow. This implies that frequencies beyond 1000 Hz also pass through. Additionally, in some figures, the low frequency component is arbitrarily represented at an even higher 2 kHz. EFP and LFP are alternately employed, and there is a lot of ambiguity with reference to notation and terminology. The specific choice for cutoff frequencies is of utmost importance because the authors' central claim about the contributions of (high frequency) action potentials propagating through the axon to LFPs. How dependent is the filter structure for the conclusions drawn by the authors? Would the conclusions be different if the authors employed, say a 20th order Butterworth filter with a 250 Hz cutoff?
We thank the reviewer for pointing out the ambiguity of our usage of the term LFP and EFP. In the revised manuscript, we now refer to the lowfrequency component as the "lowpass filtered EFP" instead of “LFP”, but point out that it is calculated in a way somewhat similar to the LFP, except for our usage of a relatively high cutoff frequency.
The reason that we chose relatively high cutoff frequencies in the kHz range was that our example simulations are meant to resemble the auditory brainstem, which is also where we obtained our data from. The optimal cutoff frequency depends on the system being studied. The auditory brainstem operates on short time scales, and the spectral frequencies are correspondingly high. In simulations with slower pulse and spike dynamics, the time scales would be slower, and frequencies would be lower.
We thus provide a general framework that can be applied to many systems and time scales, and an estimate for contributions in different scenarios given in Equation 1. Here one can see how spatial and temporal widths of the pulse, spike, and projection zone need to be modified to obtain dipolar contributions at any timescale. Furthermore, we added a paragraph to the subsection "Relevance to the interpretation of electrophysiological recordings" of the Discussion to clarify our chosen cutoff frequency.
To address the question of spectral leakage due to low filter order, we repeated the simulations for Figure 2 with a 20th order Butterworth filter (see Author response image 2), and found a very similar response to that from the 3rd order filter. This leads us to conclude that spectral leakage is not affecting the findings. In the revised manuscript, we added a note on the filtering in Materials and methods at the end of the "Multicompartment model of axons".
Finally, we now filter in the same frequency bands throughout the manuscript (low frequency < 1 kHz, high frequency > 2.5 kHz). The detailed filter structure was not important for our conclusions.
In describing the model fit, it would be valuable to describe the most important assumptions made.
In the fourth paragraph in the Results subsection (about the barn owl data), we added a description of the model assumptions when describing the fit.
Also other details of modeling need attention:
a) Location of nodes of Ranvier: In the case where there were multiple axons, were the nodes of Ranvier placed at the same spatial axis?
In the axial direction, the nodes were placed at randomized locations in all cases. They did not line up between multiple axons. This was achieved by randomizing the length of the first myelinated segment. We clarified this in the Materials and methods section.
In the simulations where the axons were placed along a line with no radial spread (all except Figure 2), the nodes were consequently also all on this line. This follows from the description in Materials and methods (section "Multicompartment model of axons"), which was improved to clarify this point.
b) "Note that in order to avoid confounding effects, the horizontal distances between axons in Figure 1CE are for illustration only; all collaterals were simulated to lie on a straight line." Isn't this an oversimplification that artificially boosts contributions of the transmembrane currents from different axons to the EFP? What happens if this is changed to reflect an anatomical organization where horizontal distances do exist?
Putting axons on a straight line has the effect of slightly underestimating the axonal contributions at all except very small (relative to the radial spread) distances from this line. There are two reasons for this:
First, a recording electrode at 100 µm from the trunk could randomly come arbitrarily close to a single fiber, leading to a very high contribution. This is not possible in our simplification. Second, if the currents in two fibers at the same depth but different radial offsets have a current with opposite polarity due to the random nature of spike generation, this would lead to an instantaneous nonzero dipole moment in the radial direction. This spurious dipole moment is different from the deterministic dipole moment discussed in the manuscript, but would lead to a stronger fluctuating contribution from the axons at increasing radial distances when compared to our simplified model (see Author response image 1). We thus consider our simplification to be a conservative one, helping us understand the axial dipole that we are interested in more clearly.
c) What was the structure of the axon at the bifurcations – was the surface area larger at bifurcations? What was the location of the terminal with reference to the nearest node?
We have clarified the corresponding section “Multicompartment model of axons” in Materials and methods. Briefly, the surface area of the bifurcation was not larger than that of the remaining axon, but was formed by a normal node segment. To avoid propagation failures at branching points, we increased the ion channel densities throughout the axon compared to those used by Simon, Carr and Shamma, 1999 in a nonbifurcating axon (see new Table 1). The nearest node to the termination was determined by the normal succession of myelinated and node segments in the axon, with the last passive segment being truncated to obtain the desired overall length of the axon.
What would be the impact of terminating axons with a bouton with appropriate dimensions containing all ion channels (calcium, sodium and potassium especially) they are supposed to express? How would the bouton geometry (Gydikov et al., 1986 and Gydikov and Trayanova, 1986) and ion channel densities there, along with changes in each parameter mentioned here, alter the conclusions?
We added a new subsection titled "Relationships to more detailed biophysical models" to the Discussion addressing this point. We argued that including all ion channels in small boutons would not alter the conclusions.
d) The authors assume that the myelin sheath does not contribute to the EFPs. However, given the several channels and receptors that these glial cells express, and given the active roles they play during action potential propagation, is it reasonable to assume that they don't contribute to the EFPs? What would be the contribution of the capacitative current that is consequent to the dielectric myelin between two parallel conducting media, with a highfrequency action potential propagating across? Is there any contribution of the 1/f^{n} scaling that is generally observed in LFPs, especially given the absence of intrinsic dendritic filtering in the analysis here (Nature Reviews Neuroscience reviews of Buzsaki et al., 2012, Einevoll et al., 2013)? Could the capacitive current (with reference to the myelinated compartments), however small, contribute to such filtering?
We addressed this point in the new subsection "Relationships to more detailed biophysical models" of the Discussion. We agree that the myelin sheath could contribute to the EFP but we think that this does not compromise our main findings. In contrast, a 1/f^{n} scaling might boost far field dipolar axonal EFPs.
e) Please provide details of how the infinite axon was simulated.
We added an explanation to the first paragraph to the Materials and methods subsection "Multicompartment model of axons" and clarified the terminology at the first usage of this term.
https://doi.org/10.7554/eLife.26106.010Article and author information
Author details
Funding
National Institute on Deafness and Other Communication Disorders (DC00436)
 Catherine Emily Carr
German Ministry for Education and Research (01GQ1001A)
 Richard Kempter
Collaborative Research in Computational Neuroscience Program (01GQ1505A)
 Richard Kempter
National Science Foundation (1516357)
 Catherine Emily Carr
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
This work was supported by the German Federal Ministry for Education and Research Grants 01GQ1001A and 01GQ0972, NIH DCD 000436 and USAmerican Collaboration in Computational Neuroscience ‘Field Potentials in the Auditory System’ as part of the NSF/NIH/ANR/BMBF/BSF Collaborative Research in Computational Neuroscience Program, 01GQ1505A. The authors wish to thank Anna Kraemer for help with animal handling and surgery; and Martina Michalikova and Tiziano D’Albis for helpful comments on a draft of this manuscript.
Ethics
Animal experimentation: This study was performed in accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. All of the animals were handled according to approved institutional animal care and use committee (IACUC) protocols R1314 and R1611 of the University of Maryland. All experiments were nonrecovery, with surgery performed under urethane anesthesia."
Reviewing Editor
 Mark CW van Rossum, University of Edinburgh, United Kingdom
Publication history
 Received: March 6, 2017
 Accepted: September 1, 2017
 Accepted Manuscript published: September 5, 2017 (version 1)
 Version of Record published: September 27, 2017 (version 2)
Copyright
© 2017, McColgan et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 989
 Page views

 188
 Downloads

 8
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.