Ephaptic coupling in white matter fibre bundles modulates axonal transmission delays

Axonal connections are widely regarded as faithful transmitters of neuronal signals with fixed delays. The reasoning behind this is that extracellular potentials caused by spikes travelling along axons are too small to have an effect on other axons. Here we devise a computational framework that allows us to study the effect of extracellular potentials generated by spike volleys in axonal fibre bundles on axonal transmission delays. We demonstrate that, although the extracellular potentials generated by single spikes are of the order of microvolts, the collective extracellular potential generated by spike volleys can reach several millivolts. As a consequence, the resulting depolarisation of the axonal membranes increases the velocity of spikes, and therefore reduces axonal delays between brain areas. Driving a neural mass model with such spike volleys, we further demonstrate that only ephaptic coupling can explain the reduction of stimulus latencies with increased stimulus intensities, as observed in many psychological experiments.


Author summary
Axonal fibre bundles that connect distant cortical areas contain millions of densely packed axons. When synchronous spike volleys travel through such fibre bundles, the extracellular potential within the bundles is perturbed. We use computer simulations to examine the magnitude and shape of this perturbation, and demonstrate that it is sufficiently strong to affect axonal transmission speeds. Since most spikes within a spike volley are positioned in an area where the extracellular potential is negative (relative to a distant reference), the resulting depolarisation of the axonal membranes accelerates the spike volley on average. This finding is in contrast to previous studies of ephaptic coupling effects between axons, where ephaptic coupling was found to slow down spike propagation. Our finding has consequences for information transmission and synchronisation between cortical areas.

Introduction
Signal processing and transmission in neuronal systems involves currents flowing across neuronal cell membranes. Due to the resistance of the extracellular medium, such transmembrane currents generate extracellular potentials (EPs), also called local field potentials (LFPs). The sources of EPs are synaptic currents, action potentials, calcium spikes and voltage-dependent intrinsic currents [1]. Neurons can therefore interact with their neighbours by changing the electric potential of the extracellular medium (and hence the membrane potential of their neighbours) without forming synapses. Such interaction is termed ephaptic interaction or ephaptic coupling [2][3][4]. Since EPs generated in the cortex are generally of the order of 100μV [5] and therefore small in comparison to neuronal threshold potentials, the influence of EPs on neural computation is often regarded as negligible. EPs can be measured with intracranial electrodes and are used as a proxy for the underlying neuronal activity [6][7][8][9]. Seminal experiments by Katz and Schmitt [10], Rosenblueth [11], Arvanitaki [2] and Marrazzi and Lorente de Nó [12] have demonstrated that action potentials travelling along parallel axons can interact with each other if the extracellular medium is highly resistive. They demonstrated that action potentials with an initial offset would resynchronise, and also slow each other down. Furthermore, action potentials could be initiated in passive axons by action potentials travelling in a nearby axon. Several studies have reproduced these effects using computational models [13][14][15][16][17][18][19][20]. However, the experimental setup is such that the axons are placed into a highly resistive medium (either paraffin oil [10], or moist air [11]) in comparison to the intracellular medium, and the computational models assume that axons are embedded within a finite-sized extracellular medium. The latter would be justified by the presence of epineuria or perineuria, which is tissue restricting the extracellular space around axons in the peripheral nervous system. Both, however, are unlikely scenarios for axonal fibre bundles within the brain: the extracellular medium is only about three times more resistive than the intracellular medium, and axons in the central nervous system are not wrapped by epineuria and perineuria that would justify the 'cables within a cable' approach. For these reasons, the amplitude of extracellular potentials around spike carrying axons should be small, and ephaptic coupling should not play a significant role between individual pairs of axons within axonal fibre bundles in the brain. However, we hypothesise that collective interaction between multiple axons affects axonal signal transmission.
We test our hypothesis by introducing a modelling framework in which EPs modulate spike thresholds, and hence spike propagation velocities. We first determine the EPs generated by action potentials in single axons, which can be computed using the axial profile of an action potential ( Fig 1A). The importance of computing the EPs lies in the fact that they perturb the membrane potential of a passive fibre (Fig 1B). This is then followed by the computation of EPs generated by spike volleys in fibre bundles (Fig 1C and 1D). As axon bundles contain millions of axons, we compute the cumulative effect of spike volleys at the macroscopic scale in axon bundles with diameters of several millimeters. The results of this analysis are used to build a spike propagation model, in which each spike travels with a velocity that is determined by structural parameters of the axon and the extracellular potential. This model is then coupled into a neural mass model to investigate in-silico the latency-intensity relationship of sensory stimuli, and the role of ephaptic coupling.
The spike propagation model that we propose in this article constitutes a strong simplification of biophysical models for spike propagation in (myelinated) axons. Rather than numerically computing the membrane potential along the axon, we asign each spike a position on an axon that changes in time according to its propagation velocity. In the absence of any external perturbations, the propagation velocity is constant along the axon and scales with its structural parameters. For example, the propagation velocity increases linearly with the axon diameter. The effect of extracellular potentials is then modelled by a coupling function that rescales the propagation velocity as a function of the EP at the spike position. The EP in the fibre bundle is computed based on core conductor theory. Each spike that travels along the fibre bundle is asigned a characteristic spatial profile, the length of which scales linearly with the propagation velocity. As this model is a strong simplification of standard biophysical models, we use a computationally feasible scenario where a synchronous spike volley interacts with itself as a test bed to calibrate the coupling term in the simplified model.
There is no direct evidence how spike propagation in axon bundles is affected by ephaptic coupling effects. We therefore present indirect evidence based on psychological experiments that investigate the relationship between the intensity of a sensory stimulus and the latency between stimulus presentation and the maximum response of the evoked potential. Such experiments have been performed for a range of sensory stimuli, including visual [21], auditory [22][23][24][25], and nociceptive stimuli [26,27]. For auditory stimuli, the first maximum (P1) is observed about fifty milliseconds after stimulus presentation, and the drop between low-intensity and high-intensity stimuli is of the order of ten milliseconds [23]. All neuronal signals, including sensory stimuli, have to pass through axon bundles to reach cortical areas. Consequently, we set up a model system in which an axon bundle is coupled to the Jansen-Rit neural mass model, which is capable to produce evoked potentials. A crucial assumption we make here is that sensory stimuli are coded as highly synchronised spike volleys, or sequences thereof. Such spike volleys could be generated by sensory neurons that encode rapid changes in sensory stimuli, or by the interaction between excitatory and inhibitory neurons. Computing the extracellular potential (EP) generated by a volley of spikes. A: An action potential, as expressed by the membrane potential V m along the axial dimension z, generates an EP that varies with z and the distance from the axon d. B: An action potential in an active axon perturbs the membrane potential of a passive axon via the EP. C: We consider spike volleys travelling along axonal fibre bundles, and D: infer from the EP the cumulative effect on the membrane potential of a passive axon. Computational studies have demonstrated that cortical circuits are capable of generating highly synchronised spike volleys with millisecond duration [28,29].

Extracellular potentials around single axons
First we computed the EPs generated by action potentials in single axons. We used the line approximation [30][31][32], given that the diameters of axons are several orders of magnitude smaller than the diameter (or the general lateral dimensions) of axonal fibre bundles: V 00 ðz 0 Þ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Here, ϕ is the EP, z is the axial dimension, d the distance from the axon, t is time, σ i and σ e are the intracellular and extracellular conductivities, a is the axon radius and V@(z) is the second derivative (curvature) of the membrane potential V. The EP is computed for different approximations of the spatial profile of an action potential, which include a piecewise linear and a piecewise quadratic approximation of spike profiles, but also spike profiles generated by a biophysical model [33] (Fig 2A-2C). The advantage of the piecewise approximation of the action potential profile is that the EP can be computed analytically (see the Methods section for details). The EP obtained from the biophysical model is computed numerically. For all the profiles we find that the maximum amplitude of the EP is of the order of microvolts (Fig 2D-2F), and at long distances d the EP decays with d −3 (Fig 2G-2I), akin to electric potentials of quadrupoles.

Extracellular potentials in fibre bundles
To compute the effect of multiple action potentials in a fibre bundle, we assumed that a completely synchronous spike volley travels through the fibre bundle. The fibre bundle was arranged as a set of concentric rings of axons, as shown in Fig 3A. The reference point to compute the EP was set at the centre of the axon bundle. We computed the EP for an increasing number of spikes, beginning with six spikes in the innermost ring of axons, then 18 spikes in the two innermost rings, and successively increasing the number of rings in which all axons carry action potentials ( Fig 3A). The maximum number of rings considered in this setup was 10 4 , which corresponds to a fibre bundle diameter of 10mm if the diameter of the uniform axons is 0.5μm. This fibre bundle contains approximately 3 × 10 8 axons, similar to the number of axons in the corpus callosum [34]. Increasing the active area (see Fig 3B for a macroscopic representation) yielded a longitudinal profile of the EP that saturated at large diameters ( Fig 3C). Interestingly, the profile is approximately proportional to −V(z), with V(z) being the spatial profile of the action potential ( Fig 3D). In the Methods section we demonstrate that this profile can be computed analytically, to a very good approximation, by the following expression: Here, g represents the ratio between the axon diameter and fibre diameter (axon plus myelin), commonly referred to as g-ratio. The relative size of the volume occupied by fibres, the fibre volume fraction, is represented by the quantity ρ, and P is the radius of the fibre bundle.
Next, we investigated how the EP changes with the position of the reference point, i.e. the point in the cross-sectional plane at which the EP is computed (Fig 4A). We found that the amplitude and longitudinal profile remained nearly unchanged, even if the reference point is close to the surface, as shown in Fig 4B and 4C. More specifically, the decrease of the amplitude is less than ten percent when the reference point is moved from the centre of the fibre bundle to 0.8 bundle radii away from the centre. Closer to the surface, the drop in amplitude is more marked. Outside of the bundle, while moving the reference point further away from the centre the EP drops rapidly, and at sufficiently large distances the drop in amplitude is proportional to d −3 . We take this as evidence that the EP at the centre of the bundle is characteristic for the EP across the entire cross-section of the fibre bundle, i.e. we assume the EP is uniform in the radial direction.
We consider spike volleys that engage all axons in the fibre bundle, which leads to EPs with amplitudes of order 100mV, as can be seen in Figs 3C and 4B. This is certainly an unphysiological scenario, since it is unlikely that all axons in a fibre bundle carry perfectly synchronised action potentials, and because such large EPs would certainly disrupt signal transmission in the participating axons. However, it is plausible that a (sufficiently synchronous) spike volley   Alternatively, one may consider a spike volley that is not perfectly synchronised, i.e. the spikes are distributed in space due to varying emission times. To illustrate the effect of such a spatial distribution, we draw spike positions randomly from a uniform distribution of varying width Δz. This spatial distribution can be associated with a temporal distribution via the relation Δz = vΔt, where v is the (intrinsic) propagation velocity of the uniform axons. In Fig 5 we show how increasing the active area affects the EP for different Δz. It can be seen that the maximum amplitude decreases with increasing Δz, and for sufficiently wide spike volleys the largest amplitude of the EP occurs near the edges of the spike volley instead of its centre.

A model for spike propagation
In addition to studying EPs generated by spike volleys in axonal fibre bundles, we are interested in the effect that EPs have on axonal signal transmission. Since the membrane potential is measured as the difference between intracellular and extracellular potential, a change of the EP implies a change of the membrane potential. For example, if the EP decreases, then the membrane potential increases, i.e. the membrane is depolarised. We assume that the EPs are not compensated by transmembrane currents, or if such processes occur, that these processes are too slow to be relevant for short spike volleys.
We begin the modelling procedure by setting up a fibre bundle with N axons, each of which has a diameter drawn from a shifted alpha distribution that was chosen to closely fit the results by Liewald et al. [35] (Fig 6A). For numerical purposes we set the number of axons N between 10 3 and 10 4 . A realistic fibre bundle contains many more axons, likely by several orders of magnitude. Conceptually, each of our model axons therefore represents a large number of axons with identical properties, but evenly distributed across the cross-sectional area. The fibre bundle is also endowed with macroscopic properties, namely the length and radius of the fibre bundle.
To test the transmission properties of a fibre bundle, we set up a spike volley with spike times drawn from a uniform distribution. The spike times define when the action potentials are generated at the proximal end of the bundle ( Fig 6B). The propagation of spikes along the axon is determined by a spike propagation model that is described in the next paragraph. The spike volley then reaches the distal end of the fibre bundle ( Fig 6C). Due to the distribution of axon diameters, this process results in a distribution of transmission delays ( Fig 6D). If the position of a spike is known, one can determine the EP generated by this spike. Since each model axon represents a large number of biological axons, we do not use the expression for single axons (Eq (1), but the one for the cumulative EPs generated by spike volleys (Eq (2)). The EP generated by a spike is thus the EP shown in Fig 3C, divided by the volume fraction occupied by the model axon. In this way, one can compute the spatial profile of the EP generated by a spike volley, see Fig 6E.
The spike propagation model tracks the position of an action potential along the fibre bundle, which is determined by the leading edge (rising phase) of the action potential. For the linear and quadratic approximations of the spike profile, the position is defined by the point where the membrane potential first deviates from resting potential. In the absence of perturbations by non-zero EPs, the velocity is constant along the fibre bundle. Therefore, the position of a spike can be tracked by multiplying the intrinsic velocity (determined by structural parameters of the axon) with the time elapsed since the spike was generated. The velocity of a spike is also determined by a putative spike threshold, which might be interpreted within a spike-diffuse-spike framework [19]. It has been demonstrated, using some simplifying assumptions, that the spike threshold V thr can be related to the activation delay Δt between two subsequent nodes of Ranvier by some nonlinear function, and therefore to the velocity of a spike [19]. In the presence of EPs, the membrane potential, and therefore the propagation velocity, is perturbed. The perturbation of the membrane potential can be intepreted as a perturbation of the spike threshold. If the membrane is depolarised (hyperpolarised) by the EP, then the spike threshold is effectively lowered (raised). For simplicity, we assume a linear relationship between V thr and Δt (Fig 6F). This results in the following relationship between the perturbed propagation velocity v and the EP: Adjusting the prefactor γ allows for the calibration of the spike propagation model with a more detailed biophysical model, which we demonstrate next.

Model calibration with biophysical model
An axon can be regarded as a core-conductor, and the spatio-temporal evolution of the membrane potential V can be described by the following cable equation: The term on the left hand side describes capacitive trans-membrane currents. The first term on the right hand side describes changes in axial currents inside the axon, and the second and third term on the right hand side describe resistive currents across the axonal membrane and the myelin sheath. The second term represents passive currents, and I(V) represents voltagegated currents, as described by the Hodgkin-Huxley model. Eq (4) describes the scenario when the extracellular potential ϕ e is zero, in which case the membrane potential V equates the intracellular potential ϕ i . If the EP is not zero, then the appropriate equation for the resulting The EP now affects the resistive currents (second and third term on the r.h.s.), as well as the capacitive current (fourth term on the r.h.s.).
We focus the calibration effort on a computationally feasible example. We consider a synchronous spike volley in a bundle composed of identical axons. The spike volley consists of spikes in one percent of all fibres. As all spikes are identical, Eq (36) is representative for all spike-carrying axons. The EP is calculated numerically at each time step from the resulting profile of the membrane potential using Eq (2). We vary the axon bundle diameter and record the change in the propagation velocity for the biophysical model and the spike propagation model (Fig 7).
To fit the spike propagation model to the biophysical model, we adjust the coupling parameter γ and the standardised spike profile. Specifically, we adjust the length of the spike profile and the position of the maximum. The parameter γ determines the amount of relative slowing down that is observed in both models, whereas the profile of the model spikes determines how this effect changes with varying bundle diameters (Fig 7).

Effect of extracellular potentials on transmission delays
The spike propagation model allows us to test the consequences of ephaptic coupling via EPs in a macroscopic fibre bundle. We investigate the dynamics of spike volleys with and without ephaptic coupling, and the resulting differences in axonal delays. There are several structural parameters that we keep fixed for simplicity, such as the fibre volume fraction (80% [36]), the fibre length (10cm), and the distribution of axon diameters. The spikes are generated at the proximal end of the fibre bundle, with spike times drawn from a uniform distribution. The width of this distribution determines the duration of a stimulus, and the number of spikes determines its intensity. We record the arrival time when a spike reaches the distal end of the axon, and the difference between the arrival time and the time the spike was initiated at the proximal end constitutes the axonal delay.
We first investigated how axonal delays are affected by ephaptic coupling, and focused on the mean of the delay distribution. In the presence of ephaptic coupling, we observe a decrease of the mean axonal delays as the stimulus intensity is increased (solid lines in Fig 8). In the absence of such coupling, the mean axonal delays remain constant (dashed lines in Fig 8). The stimulus duration is set to either 1ms, 2ms and 3ms, and the bundle diameters are varied between 2mm and 8mm. The mean axonal delays drop nonlinearly with increasing intensity in the presence of ephaptic coupling, but remain unchanged in its absence (Fig 8A-8C). At full intensity (100%) and with ephaptic coupling, the mean axonal delays drop from 35ms to 20ms as the diameter of the fibre bundle is increased from 2mm to 8mm if the stimulus duration is 1ms (Fig 8A). At 2ms stimulus duration, the mean axonal delays drop from 36ms (unchanged) to 28ms with increasing diameter of the fibre bundle (Fig 8B). In other words, the mean axonal delay decreased by up to 40% at full stimulus intensity.
Next, we explored how the standard deviation of axonal delays (a measure for its dispersion) behaved in the presence of ephaptic coupling. We found that its qualitative behaviour is PLOS COMPUTATIONAL BIOLOGY different from the mean of the axonal delays (Fig 8D-8F), with an initial decrease and a subsequent increase in the standard deviation.
Finally, we incorporated the axon bundle into the Jansen-Rit neural mass model [37]. The arrival of each spike at the distal end generates a current that is injected into the neural mass model. The response latency is determined by the time difference between stimulus onset and the maximum response of the neural mass model. This results in increased latencies as the stimulus duration is increased. However, in the presence of ephaptic coupling, we observe again a decrease in latencies as the stimulus intensity is increased, whereas in the absence of ephaptic coupling the decrease is only marginal (Fig 8G-8I). Regardless of stimulus duration, at full stimulus intensity ephaptic coupling reduces the response latency by up to 8ms, which corresponds to a reduction by approximately 15%.

Discussion
The key finding of our study is that spike volleys generate EPs with sufficiently large amplitudes to modulate axonal delays. Specifically, the mean delay of a spike volley decreases as the number of spikes in the spike volley is increased. Therefore, our results suggest that varying the amplitude of a neuronal signal can adjust its delay. Using a neural mass model, we have demonstrated that the decrease of axonal delays translates into the decrease of response latencies as the stimulus intensity is increased.
We have calibrated the spike propagation model using a biophysical model, by comparing the velocity change of a spike volley within a fibre bundle composed of identical axons. Here we observed the opposite effect: In the presence of ephaptic coupling, the spike volley slowed down. This is in line with previous numerical studies which investigated ephaptic coupling effects between a small number of identical axons. There, ephaptic coupling led to a synchronisation of the spikes within a volley, and a concurrent deceleration of the spikes. The acceleration of spike volleys that we observe in fibre bundles with distributed axon diameters can be attributed to the dispersive effect, which results from the axon diameter distribution and prevents synchronisation of a spike volley. As we show in Fig 5, an axially distributed spike volley causes primarily depolarisation within the fibre bundle, which then results in the acceleration of the majority of spikes within the volley. Therefore, our results do not contradict previous studies, but generalise previous modelling approaches. Nevertheless, the spike propagation model that we devised here required several assumptions that we are going to discuss in more detail.
We computed the EP using the line approximation (i.e. the axon is assumed to be infinitely thin), which has been demonstrated to be very accurate [30]. We further assumed that the axon bundle is large, circular, homogeneous, and densely populated with axons. The latter is justified by electron micrography studies which suggest that only a small fraction of an axon bundle is made up of extracellular space [35,38]. Since axonal membranes and the myelin sheaths have a much larger resistivity compared to the extracellular medium, electric currents can only pass through the extracellular medium. We assumed that the medium is homogeneous and that the effective conductivity of the fibre bundle is the conductivity of the extracellular fluid multiplied by the relative size of the extracellular space. This calls for more detailed simulations of the spread of EPs with spatial heterogeneity taken into account. For mathematical convenience, we chose the fibre bundles to be large with circular cross sections. Realistic fibre bundles are indeed large, but often show a more sheet-like morphology [39]. An open question is whether this morphology influences the effect of EPs within our framework (a recent study used coupled axons with FitzHugh-Nagumo dynamics to demonstrate ephaptic coupling effects in sheet-like bundles [20]). Furthermore, we ignored possible effects due to the axonal microstructure. We assumed the axonal membrane to be smooth (effectively a homogenised axon [40]), and that therefore nodes of Ranvier are not relevant as point sources. This is certainly the case at large distances from the axon, as can be seen in Fig 2F. However, at close proximity such effects would be relevant, as the EP at a node of Ranvier can reach several hundreds of microvolts. It is unknown whether nodes of neighbouring axons are sufficiently aligned to affect action potential generation in such a manner. As oligodendrocytes can myelinate multiple axons [41,42], it is conceivable that neighbouring axons show some degree of alignments, in which case it would be possible to observe ephaptic coupling effects in much smaller fibre bundles, provided the spike volleys are sufficiently synchronised. Furthermore, we have demonstrated that computing the EP in a fibre bundle can be very well approximated by Eq (2). This expression depends only on the membrane potential V(z) and a convolution thereof. As can be seen from Fig 2, the spatial profile of the membrane potential is quite smooth along the axon, which is due to the fact that the length of a characteristic spike is by at least one order of magnitude larger than the length of a myelinated segment.
To demonstrate the effect of EPs on axonal delays we used a strongly simplified model for spike propagation. This model assumes that the spike velocity resulting from the axon morphology is known, and that this velocity is perturbed by the EP. It does not contain possible compensation effects arising from Hodgkin-Huxley dynamics (i.e. subthreshold currents that repolarise the axonal membrane), and further studies are required using the Hodgkin-Huxley framework to confirm our results. While we have used a simple scenario to calibrate the spike propagation model, to reproduce all of our findings within the Hodgkin-Huxley framework would be computationally extremely expensive, and is therefore beyond the scope of the present study.
We have incorporated the axon bundle model into the Jansen-Rit neural mass model to build a model system for primary sensory information processing, and to investigate the relationship between stimulus intensity and response latency. Psychological experiments across different sensory modalities yield the same qualitative relationship, whereby the latency decreases with increasing intensity [21][22][23][24][25][26][27]. Such experiments typically measure the delay between stimulus presentation and the first maximum or minimum of the neural response measured electrographically. To replicate this experimental design, we measured the time difference between stimulus onset, i.e. the start of the spike volley, to the first maximum in the response of the Jansen-Rit model. Interestingly, only the presence of ephaptic coupling could explain the latency-intensity relationship. We are aware that the Jansen-Rit model is a fairly simple representation of a cortical microcircuit, and that other nonlinear processes not taken into account in its derivation may also reduce the response latency with increased stimulus intensity, such as oscillation-mediated information transmission [43]. Nevertheless, our modelling approach suggests that ephaptic coupling effects play a role in neural responses.
While there is such implicit evidence, further experimental studies are necessary to test our hypothesis. The experimental design would be highly invasive, since the EPs drop rapidly with distance outside fibre bundles. Animal experiments have already demonstrated the possibility to record EPs within axonal fibre bundles [32,44]. An interesting test bed could also be a delay analysis within stimulation-response paradigms used in epileptic patients to determine the seizure focus [45,46].
If such activity-dependent (or rate-dependent) delays occur in fibre bundles, then one may speculate as to their putative role in information processing. Since axonal delays are in general quite small (about 30ms in a 10cm long fibre bundle), the main effect should be on fast oscillations. It is indeed tempting to propose that such variable delays may have an effect on longrange gamma synchronisation, and that synchronisation patterns can be flexibly switched by changes in the amplitude of the transmitted spike volleys. We have found that ephaptic coupling can decrease delays by up to 10ms, which would be one half of the period of a gamma cycle at 50Hz. It has been demonstrated that delays are critical in shaping the functional architecture of the brain [47][48][49], and ephaptic modulation of such delays could therefore play a role in flexibly synchronising distant brain areas.

Methods
In this section, we describe the mathematical framework underlying our study. We use a detailed microscopic description of the interaction between axonal fibres, and a leading-edge approximation which reduces the computational effort, but retains key properties of the interaction. We investigate a fibre bundle in which axons are coupled by the EPs generated by spikes. We first show how to compute the EPs generated by single spikes and spike volleys, and then present the framework of the spike propagation model.

Single fibre
In an open fibre bundle, the EP is determined by currents entering and leaving an axon. The radial currents around an axon can be inferred from the spatial profile of an action potential [32] given by the cable equation: with I(z, t) being the radial (trans-membrane) currents, V(z, t) the membrane potential, a the axon radius and σ i the conductivity of the intracellular medium. The axial dimension is represented by z, and time by t. The EP, denoted by ϕ for the single spike, can then be computed from the radial current via: Iðz 0 ; tÞ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi with σ e being the extracellular conductivity, and d the radial coordinate measuring the distance from the axon. Inserting Eq (6) into Eq (7)  V 00 ðz 0 Þ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi This integral is the convolution of the curvature of the action potential profile with the kernel G(z − z 0 ) = ((z − z 0 ) 2 + d 2 ) −1/2 . In general, this integral has to be evaluated numerically. In order to obtain an analytical solution, we approximate the shape of an action potential by piecewise linear or piecewise quadratic functions. Piecewise approximation of action potential profile. In general the profile of an action potential has to be determined either numerically, or using spike-diffuse-spike formalisms. In the former case it is impossible to parameterise the profile, and in the latter the analytical expressions are still prohibitive to follow through with the calculations of the EP. Therefore, we present a formalism which approximates the profile of an action potential with either piecewise linear or piecewise quadratic functions. This method can be extended to arbitrary polynomial expressions, and is similar to curve-fitting with splines.
piecewise linear approximation The simplest approximation of an action potential is given by two linear functions on two consecutive intervals, describing the rising and the falling phase of the action potential, respectively: The first derivative of this approximation is piecewise constant with discontinuities at z = z 0 , PLOS COMPUTATIONAL BIOLOGY z = z 1 and z = z 2 . The second derivative is therefore It is then straightforward to compute the EP: �ðz; d; tÞ ¼ s i a 2 4s e V max z 1 À z 0 1 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi � 1 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi piecewise quadratic approximation For the piecewise quadratic approximation, we divide the AP profile into three segments: Given z 0 , z 1 , z 2 , z 3 and V max , there are four unknowns a 1 , a 2 , a 3 and z max . To ensure a smooth profile, we impose boundary conditions that assume V(z) is smoothly differentiable, i.e. Vðz ! z þ After some manipulation, we obtain: The second derivative of the spatial profile is piecewise constant: The EP is then found to be �ðz; d; tÞ ¼ s i a 2 2s e a 1 ln ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi d 2 þ ðz À z 1 Þ 2 q þ z 1 À z ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi À a 2 ln ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi d 2 þ ðz À z 2 Þ 2 q þ z 2 À z ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi þa 3 ln ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi d 2 þ ðz À z 3 Þ 2 q þ z 3 À z ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi

Fibre bundle
Bundle with identical axons. To compute an upper boundary of the EP produced by multiple action potentials, we assume that a perfectly synchronous spike-volley travels through a dense fibre bundle. All axons in this fibre bundle have the same diameter and are arranged in concentric rings (Fig 3A). At the centre of an empty grid position we compute the EP by summing ϕ at distance (2n + 1)a of 6n axons, with n ranging from 1 to N, with N large: Although ϕ is approximately 20μV at the surface of an isolated axon, in a fibre bundle the combined effect can lead to EPs of many mV. Interestingly, we find that for large enough fibre bundle diameters the profile of the cumulative EP is almost proportional to the profile of the generating action potentials. We give a mathematical explanation for this next. Analytical solution. The cumulative EP at the core of an axon bundle is computed with the following integral, with P being the axon bundle diameter, and O(a) being the cross-sectional area occupied by an axon with radius a. We set O(a) = πa 2 /(ρg 2 ), where g is the g-ratio and ρ is the fibre volume fraction. Inserting Eq (1) into Eq (23), and solving the integral over ρ, results in V 00 ðz 0 Þ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Integration by parts then yields EPðz; PÞ ¼ s i g 2 r 2s e Z 1 À 1 V 0 ðz 0 Þ ðz À z 0 Þ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ðz À z 0 Þ 2 þ P 2 q À sgnðz À z 0 Þ Next, we use the approximation ðz À z 0 Þ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ðz À z 0 Þ 2 þ P 2 q � sgnðz À z 0 Þ 1 À e À jzÀ z 0 j=P À � ; ð26Þ which leads to Using integration by parts, this integral ultimately yields EPðz; PÞ � À s i g 2 r s e VðzÞ þ s i g 2 r 2s e 1 P We may regard this result as the far-field approximation of EPs in axonal fibre bundles.
We note here also that in the limit P ! 0, exp(|z − z 0 |/P)/P ! 2δ(z − z 0 ), which yields EPðz; 0Þ ¼ À s i g 2 r s e VðzÞ þ s i g 2 r s e VðzÞ ¼ 0: Eq (24) suggests that the far-field approximation of the cumulative EP is independent of the axon morphology. At this point, however, we have not taken into account that axons of different diameters transmit action potentials at different velocities, and that therefore the spatial profile widens with increasing action potential velocity.

Spike propagation model
The two most common ways to model axonal signal transmission are either Hodgkin-Huxley type dynamics embedded in a core-conductor model, or simpler spike-diffuse-spike approaches. Both ways allow one to determine the spike velocity as a function of electrophysiological and structural parameters. Here, we employ a much simpler model that describes the position z i of an action potential (more precisely, its leading edge or rising phase) travelling along the i th axon by one simple equation: where v i (z, t) is the velocity of the action potential as function of the axial direction z and time t. If the axon is homogeneous and does not experience spatial or temporal perturbations, then the velocity can be expressed by v i (z, t) = v i,0 , which is the intrinsic velocity of the axon, determined by its morphological and electrophysiological properties. We assume here that this velocity is known for each axon. In the absence of perturbations, one can therefore express the axonal delays by τ i = L/v i,0 , with L being the length of the fibre bundle. We set here v i,0 = αd i , with d i being the diameter of the i th axon, and α = 5ms −1 /μm. Changes in the EP lead to perturbations of the membrane potential of an axon. A negative (positive) EP effectively depolarises (hyperpolarises) the axonal membrane, and therefore increases (decreases) the propagation velocity. A convenient formalism to incorporate such changes is the spike-diffuse-spike framework, in which the spiking threshold is a parameter to explicitly describe the onset of an action potential [19]. Such thresholds can also be determined within the Hodgkin-Huxley formalism, albeit these thresholds vary with the depolarisation rate [50]. The EP can therefore be regarded as a perturbation of such a threshold. Within the spike-diffuse-spike framework, the following relationship between the spike threshold V thr and the delay of spike generation Δt between two consecutive nodes of Ranvier can be derived: see Fig 6F for a visual representation. The function f(Δt) depends on structural and electrophysiological parameters of the axon. The EP can be incorporated into the spiking threshold, V thr (z, t) = V thr,0 + EP(z, t), with V thr,0 being the uniform spiking threshold of the unperturbed axon. Via Eq (31) one can relate V thr,0 to Δt 0 of the unperturbed axon, and to its intrinsic velocity via v i,0 = l/Δt 0 , where l is the distance between two consecutive nodes of Ranvier. To further simplify our scheme we linearise Eq (31) around V thr,0 and Δt 0 : Here we made use of v = Δz/Δt(z, t), which results in Δt(z, t)/Δt 0 = v 0 /v(z, t). The parameter c denotes the relative steepness of f(Δt) around Δt 0 . For simplicity, we lump c and V thr into one parameter γ = 1/(cV thr ), which results in Eq (3). This parameter can be used as a tuning parameter for the inverse strength of the ephaptic coupling, and in the numerical simulations we set γ = 0 to represent the absence of any ephaptic interaction. The computation of the EP generated by a spike requires knowledge about its axial profile, which can be inferred from the temporal evolution of a spike at any given point along the axon. In order to obtain a computationally efficient framework, we focus on the piecewise linear approximation. The temporal profile of a spike in linear approximation is given by where t 0 is the time of onset of the spike, t 1 is the time it reaches its maximum depolarisation, and t 2 is the time when the membrane is fully repolarised. The variable v(z, t) represents the velocity of the leading edge of the spike, which can vary rapidly. We introduce the effective velocity to estimate the velocity of the entire spike, or of its centre of mass. This effective velocity is, in essence, a time-averaged quantity. We approximate the effective velocity of a spike by with τ = 1ms. The temporal profile can then be translated into the spatial profile by setting z = v eff (z) × t, and z m = v eff (z) × t m . The parameters for the spike propagation model are listed in Table 1.

Biophysical model
An axon can be regarded as a core-conductor with intracellular potential ϕ i (z, t), and extracellular potential ϕ e (z, t). The difference between in intracellular and extracellular potential results in the membrane potential V(z, t) = ϕ i (z, t) − ϕ e (z, t). The extracellular potential is subject to the model setup. Here we consider an axon bundle composed of identical axons, a fraction of which carries spikes synchronously. Because the extracellular potential is considered homogeneous within the fibre bundle (Eq 2), the dynamics of all spikes is identical, and it suffices to solve one cable equation to describe the propagation of all spikes: Since we consider myelinated axons, the radial capacitance C m and the radial resistance R m vary with the axial coordinate z, depending on whether a segment is myelinated (internode) or unmyelinated (node of Ranvier): R my else: The model axon is myelinated periodically with nodes of length l, capacitance C n , and resistivity R n , and myelinated internodes with length L, capacitance C my and resistivity R my . The intracellular resistance R ax remains constant along the axon. The term I(V) represents voltage-gated currents modelled with the Hodgkin-Huxley model. These currents only occur in nodal segments: IðVÞ ¼ g Na;f m 3 hðe Na À VÞ þ g Na;p p 3 ðe Na À VÞ þ g K;s sðe K À VÞ: The gating variables m, h, p, s obey the following dynamics: _ n ¼ a n ðVÞð1 À nÞ À b n ðVÞn; ð39Þ with n = m, h, p, s. The variables α n and β n are defined as follows: a m ðVÞ ¼ 1:86ðV þ 25:4mVÞ=ð1 À expðÀ ðV þ 25:4mVÞ=10:3mVÞÞ=ms; b m ðVÞ ¼ 0:086ðÀ ðV þ 29:7mVÞÞ=ð1 À expðÀ ðV þ 29:7mVÞ=9:16mVÞÞ=ms; a h ðVÞ ¼ 0:0336ðÀ ðV þ 118mVÞÞ=ð1 À expðÀ ðV þ 118mVÞ=11mVÞÞ=ms; The parameters for the biophysical model are given in Table 2.

Jansen-Rit microcircuit
The spike volleys represent cortical, subcortical or sensory information being transitted by axon bundles. To describe the response of neuronal circuits, e.g. cortical microcircuits, we use the Jansen-Rit model [37,51,52] and record the maximum response in the membrane potential of its pyramidal cell population. The Jansen-Rit model is composed of six differential equations: _ y 0 ¼ y 3 _ y 3 ¼ Aas½y 1 À y 2 � À 2ay 3 À a 2 y 0 _ y 1 ¼ y 4 _ y 4 ¼ AapðtÞ þ C 2 s½C 1 y 0 � À 2ay 4 À a 2 y 1 _ y 2 ¼ y 5 _ y 5 ¼ BbC 4 s½C 3 y 0 � À 2by 5 À b 2 y 2 : The diameter was adjusted to produce a propagation velocity of 4m/s. The parameters related to the Hodgkin-Huxley model are taken from Ref [33]. The values related to the variables α n , β n are not listed here, but stated explicitly in the main text. https://doi.org/10.1371/journal.pcbi.1007858.t002 Here, y 0 is the postsynaptic potential (PSP) generated by the output of the pyramidal cells at the two interneuron types, and y 1 and y 2 are the excitatory and inhibitory PSPs generated at the pyramidal cells by external stimuli p(t) and the firing activity of the interneurons. Furthermore, y 3 , y 4 and y 5 are auxiliary variables in the synaptic conversion of firing rates into PSPs, with a and b being the inverse time constants of excitatory and inhibitory synapses, and A and B the respective maximum amplitudes of the synaptic response. The PSPs are converted into firing rates by the sigmoidal function s½v� ¼ e 0 1 þ e rðv 0 À vÞ ; ð41Þ with e 0 being the maximum firing rate, v 0 the membrane potential at half of the maximum firing rate, and r sets the steepness of the sigmoid. The interaction between the different neuron types is scaled by the connectivity constants C 1 to C 4 . The parameters are chosen as in [37] (see Table 3). The external firing rate p(t) is generated by the incoming spikes, with t n being the arrival time of the n th spike, N the total number of spikes, and P = 0.1s −1 sets the synaptic coupling strength for the spike train. The membrane potential y of the pyramidal cell population is determined by the difference between excitatory and inhibitory PSPs, i.e. y = y 1 − y 2 . The response latency is then calculated as the position of the absolute maximum of y.