Extracellular potentials around single axons

First we computed the EPs generated by action potentials in single axons. We used the line approximation [30–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:

Fig 2

Spatial profiles of action potentials and their EPs.

Shown are A: the piecewise linear profile, B: the piecewise quadratic profile, and C: the profile of an action potential generated with the biophysical model. D-F: EPs corresponding to action potential profiles in A-C. G-I: Log-log plots of the EPs (absolute values) at z = 0. Black lines indicate decay with d⁻³. (The notch at d ≈ 0.3mm is due to a change of sign).

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⁴, 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⁸ axons, similar to the number of axons in the corpus callosum [34].

Fig 3

EP at the centre of a circular axon bundle due to concentric spike volleys.

A: Microscopic cross-section of a fibre bundle, with spike-carrying axons marked in blue. B: Macroscopic extension of (a), with the active area (i.e. where axons carry spikes) marked in blue. C: Waveform of a spike (top), and the resulting spatial (axial) profile of the EP at the centre of the fibre bundle. D: Cross-sections of C.

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:

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⁻³. 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.

Fig 4

EP in fibre bundle with synchronous spike volley, subject to position of reference point.

A: The reference point is moved from the centre of the fibre bundle to a position outside of the fibre bundle. B: Waveform of a spike (top), and the resulting EP plotted against the longitudinal coordinate z and the distance of the reference point from the centre. C: Cross-sections of B.

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 might engage one percent of the axons in the fibre bundle, in which case the amplitude of the EP would be of the order of 1mV.

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.

Fig 5

Increasing the length of a spike volley attenuates the amplitude of an EP.

The EP is shown for varying bundle diameters and z. We steadily increase the width Δz of the spike volley from A: Δz = 0mm, to F: Δz = 50mm.

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³ and 10⁴. 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.

Fig 6

Illustration of properties of the computational model.

A: Distribution of axon diameters sampled from a shifted alpha distribution to match experimental data [35]. B: Rastergram of spike volley generated at proximal end of fibre bundle. C: Rastergram of spike volley reaching the distal end of the fibre bundle. D: Distribution of delay times. E: Snapshot of the longitudinal profile of EP generated by a spike volley. F: The EP modulates the spiking threshold (V_thr) and therefore the delay Δt of action potential generation between two reference points (e.g. two consecutive nodes of Ranvier).

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:

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:

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

Fig 7

Comparison of the spike propagation model with a biophysical model.

A synchronous spike volley slows down as a result of ephaptic coupling in fibre bundles with identical axons. The relative change of the propagation velocity varies with the bundle diameter.

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.

Fig 8

Increasing the stimulus intensity, i.e. the number of spikes in a volley, decreases axonal transmission times and the latency of stimulus response.

A-C: Mean axonal delay with ephaptic coupling (solid) and without ephaptic coupling (dashed) for A: 1ms, B: 2ms, and C: 3ms stimulus duration. D-F: Standard deviation from the mean of axonal delay with ephaptic coupling (solid) and without ephaptic coupling (dashed) for D: 1ms, E: 2ms, and F: 3ms stimulus duration. Mean and standard deviation are computed from the distribution of delay times (cf. Fig 6D). G-I: Latency from stimulus onset to first maximum in neural mass model at G: 1ms, H: 2ms, and I: 3ms stimulus duration. Lines (shaded areas) indicate mean (1σ confidence interval) across 5 simulations. Colours indicate different bundle diameters.

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

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:

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₀, z = z₁ and z = z₂. The second derivative is therefore

It is then straightforward to compute the EP:

piecewise quadratic approximation

For the piecewise quadratic approximation, we divide the AP profile into three segments:

Given z₀, z₁, z₂, z₃ and V_max, there are four unknowns a₁, a₂, a₃ and z_max. To ensure a smooth profile, we impose boundary conditions that assume V(z) is smoothly differentiable, i.e. $$ V(z\to z_1^{+})=V(z\to z_1^{-})$$, $$ V^{\prime }(z\to z_1^{+})=V^{\prime }(z\to z_1^{-})$$, $$ V(z\to z_2^{+})=V(z\to z_2^{-})$$, and $$ V^{\prime }(z\to z_2^{+})=V^{\prime }(z\to z_2^{-})$$. After some manipulation, we obtain:

The second derivative of the spatial profile is piecewise constant:

The EP is then found to be

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 Ω(a) being the cross-sectional area occupied by an axon with radius a. We set Ω(a) = πa²/(ρg²), 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

Integration by parts then yields

Next, we use the approximation

which leads to

Using integration by parts, this integral ultimately yields

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′|/P)/P → 2δ(z − z′), which yields

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:

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:

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

Table 1

List of parameters used for the spike propagation model.

Parameter	standard values
N	104
L	10cm
α	5ms−1/μm
γ	1/180mV
Vmax	100mV
t1	0.3ms
t2	2ms
τ	1ms
ρ	0.8
σi/σe	3/(1 − ρ)
g	0.8

Unless explicitly stated, we use the parameters presented in this list. For most figures we use the standard parameters.

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:

The term I(V) represents voltage-gated currents modelled with the Hodgkin-Huxley model. These currents only occur in nodal segments:

The parameters for the biophysical model are given in Table 2.

Table 2

List of parameters used for the biophysical model.

Parameter	standard values
d	0.7706μm
L	200d
l	2μm
g	0.6
Cmy	3.6 × 10−6/ln(1/g)μFcm−2
Cn	1μFcm−2
Gmy	7.7 × 10−6/ln(1/g)mScm−2
Gn	80mScm−2
Gax	14.3mScm−1
ρ	0.8
σi/σe	3/(1 − ρ)
gNa,f	3000mScm−2
gNa,p	5mScm−2
gK,s	80mScm−2
eNa	132mV
eK	−2mV

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.

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:

Table 3

List of parameters used for the Jansen-Rit model.

Parameter	standard values
A	3.25mV
B	22mV
a	100s−1
b	50s−1
v0	6mV
e0	5s−1
r	0.56mV−1
C1	135
C2	0.8
C3	33.75
C4	33.75
P	0.1s−1

Parameters are taken from Ref [51] and used throughout the manuscript. P is adjusted such that even at low intensities a response with clear maximum is observable.

The external firing rate p(t) is generated by the incoming spikes,

The membrane potential y of the pyramidal cell population is determined by the difference between excitatory and inhibitory PSPs, i.e. y = y₁ − y₂. The response latency is then calculated as the position of the absolute maximum of y.