Modelling the effect of ephaptic coupling on spike propagation in peripheral nerve fibres

Experimental and theoretical studies have shown that ephaptic coupling leads to the synchronisation and slowing down of spikes propagating along the axons within peripheral nerve bundles. However, the main focus thus far has been on a small number of identical axons, whereas realistic peripheral nerve bundles contain numerous axons with different diameters. Here, we present a computationally efficient spike propagation model, which captures the essential features of propagating spikes and their ephaptic interaction, and facilitates the theoretical investigation of spike volleys in large, heterogeneous fibre bundles. We first lay out the theoretical basis to describe how the spike in an active axon changes the membrane potential of a passive axon. These insights are then incorporated into the spike propagation model, which is calibrated with a biophysically realistic model based on Hodgkin–Huxley dynamics. The fully calibrated model is then applied to fibre bundles with a large number of axons and different types of axon diameter distributions. One key insight of this study is that the heterogeneity of the axonal diameters has a dispersive effect, and that a higher level of heterogeneity requires stronger ephaptic coupling to achieve full synchronisation between spikes.


Introduction
Signal transmission in peripheral nerve fibres is based on the propagation of action potentials, or spikes, along the axonal membrane, which alters the electrophysiological properties of the extracellular medium and therefore influences nearby axons via so-called ephaptic coupling (Anastassiou et al. 2011;Buzsáki et al. 2012;Anastassiou and Koch 2015 Institute of Biomedical Engineering and Informatics, University of Technology Ilmenau, Gustav-Kirchhoff Str. 2, 98693 Ilmenau, Germany of a highly resistive extracellular medium, action potentials travelling in two parallel, closely spaced axons synchronise and travel at a lower velocity than in the isolated case (Katz and Schmitt 1940;Rosenblueth 1941;Arvanitaki 1942;Marrazzi and Lorente de Nó 1944). This has been reproduced in theoretical work using numerical or analytical tools (Clark and Plonsey 1970;Bell 1981;Eilbeck et al. 1981; Barr and Plonsey 1992;Binczak et al. 2001;Bokil et al. 2001;Reutskiy et al. 2003;Maïna et al. 2015;Shneider and Pekker 2015;Goldwyn and Rinzel 2016; Knösche 2019; Sheheitli and Jirsa 2020; Capllonch-Juan and Sepulveda 2020). Both theoretical and experimental work, however, have been restricted thus far to a small number of identical axons due to the experimental or computational effort (with the exception of Capllonch-Juan and Sepulveda (2020)). In contrast, peripheral nerve bundles are composed of a relatively large number of nerve fibres, which follow a wide distribution of axonal diameters (Sanders 1947;Assaf et al. 2008;Ikeda and Oka 2012;Eichel et al. 2020). As there is an approximately linear relationship between the axon diameter and the velocity of a spike (Hursh 1939;Goldman and Albus 1968;Schmidt and Knösche 2019), one can expect that structural heterogeneity counteracts the synchronisation of spikes, akin to phase-coupled oscillators, where wider distributions of the angular frequencies require stronger coupling to achieve synchronisation (Kuramoto 1984;Pikovsky et al. 2001). The main difference between phase-coupled oscillators and ephaptically coupled spikes is that in the latter case the coupling is restricted to axonal segments close to the spikes, and that the lifetime of a spike is restricted to the time interval between spike generation and the spike reaching the distal end of the axon. Therefore, spike synchronisation depends crucially on the initial conditions, as demonstrated recently for white-matter fibre bundles (Schmidt et al. 2021). For this reason, we focus our analysis on initially synchronous spike volleys, in line with previous studies (Binczak et al. 2001;Reutskiy et al. 2003).
Our aim is to systematically study ephaptic coupling effects in peripheral nerves with large axon counts and distributed morphology. For this, we devise a simplified spike propagation model (SPM) that reduces the computational effort significantly as compared to solving models based on nonlinear partial differential equations. The SPM is inspired by phase-coupled oscillators in the sense that each spike is endowed with an intrinsic propagation velocity based on the morphology of the associated axon. Furthermore, each spike can be ascribed a coupling function that characterises the spatiotemporal profile of ephaptic coupling with spikes in nearby axons. The collective ephaptic coupling generated by all spikes in a fibre bundle then perturbs the intrinsic velocity of a single spike, which decelerates if the axonal membrane is hyperpolarised, and vice versa. In the resulting numerical scheme, the position of each spike is propagated based on its perturbed velocity. In turn, the positions and velocities of the spikes determine their mutual ephaptic coupling effects.
We begin our investigation by considering an exemplary bundle containing two axons, for which we develop the theoretical basis that can be extended to larger fibre bundles. First, we compute the effect that the spike in an active axon has on a nearby passive axon, which we do analytically under the assumption that the axons are homogeneous and that the perturbation of the passive axon can be captured by the linear terms of the cable equation. As most axons in peripheral nerves are myelinated, the assumption of homogeneity implies that the specific location of nodes of Ranvier is irrelevant. We calibrate the proposed SPM using a biophysically realistic model based on Hodgkin-Huxley dynamics.
The calibrated SPM is then used to study large axon bundles systematically, especially the interplay between axon diameter distribution and the strength of ephaptic coupling, with the latter having a synchronising effect, and the former a dispersive effect. The coupling strength is determined by the fibre density (i.e. relative volume occupied by nerve fibres), and the resistivity of the extracellular medium. The aim here is to identify conditions under which synchroni-sation can occur given different axon diameter distributions and structural parameters, especially fibre density.

Ephaptic coupling between nerve fibres
In this section, we lay the theoretical basis of the SPM by considering fibre bundles with a small number of axons. First, we derive an analytical expression for the extracellular potential generated by a spike in an active axon. This is then used to compute the perturbation of the membrane potential of a passive axon (an axon that does not carry a spike) generated by a spike in a nearby active axon. The perturbation can be regarded as the ephaptic coupling between two nerve fibres. The resulting mathematical framework is then extended to the interaction between multiple axons. To reduce the computational load, we introduce a piecewise quadratic approximation of the spatial spike profile, which allows us to compute the ephaptic coupling function analytically. Lastly, we introduce the SPM and calibrate it by fitting free parameters to a biophysically realistic model based on the cable equation and Hodgkin-Huxley dynamics.

Perturbation of the extracellular medium by a spike in a single axon
Axons may be regarded as core conductors (Rall 1977;Trayanova et al. 1990;Holt and Koch 1999), whereby the formation and propagation of spikes can be described by the one-dimensional cable equation: The derivation of the cable equation is based on Kirchhoff's first law, which states that the currents are balanced at any given location inside or outside the axon. Here, the membrane potential V is defined as the difference between intracellular and extracellular potential, V = φ i − φ e . If the extracellular medium is grounded, then V = φ i and the traditional cable equation is recovered. In general, however, that assumption does not hold because of the finite size of the extracellular medium, and its finite conductivity (Tveito et al. 2017). The term on the left-hand side of Eq. (1) describes the capacitive currents across the axonal membrane and the myelin sheath, with C m being the capacitance of the membrane and myelin. The terms on the right-hand side of Eq. (1) describe resistive currents, of which the first term is the longitudinal intra-axonal resistance, which depends on the intracellular conductivity σ ax and the cross-sectional area of the axon A ax = πr 2 ax via g ax = σ ax A ax . The second term is the passive resistive current across the membrane and myelin, with g m being the radial resistance of the myelin and the mem-brane, and the third term represents voltage-gated currents described by the Hodgkin-Huxley model.
Under the assumption that the boundary of the fibre bundle is perfectly insulating, the current balance outside an axon can be formulated in the equivalent manner: Here, the (longitudinal) extracellular resistance g ex depends on the extracellular conductivity σ ex and the cross-sectional area of the extracellular space, A ex . This relationship, however, is no longer applicable if multiple axons are in the same bundle. In this case, the extracellular potential depends on the capacitive and resistive currents generated by all axons, and we therefore find the following relationship for the current balance in the extracellular medium: with N being the total number of axons. It is easy to see from Eq.
(3) that if the extracellular conductivity or the crosssectional area of the extracellular medium increases (thus increasing g ex ), then φ e decreases (which has been demonstrated numerically in Tveito et al. (2017)).
Typically, spikes are measured as the spatiotemporal profile of the depolarisation of the membrane potential V . Here, we are interested in recovering a relationship between the extracellular potential φ e and the membrane potential V . Combining Eq. (3) and Eq. (1), and using φ i,n = V n + φ e , we obtain φ e = − n g ax,n V n g ex + n g ax,n .
The terms V n and φ e denote the second spatial derivative in the axial direction of V n and φ e . Although Eq. (4) is only an implicit representation of the extracellular potential via its second spatial derivative, we show that this expression is sufficient to compute the effect of a spike onto the membrane of a nearby axon. It is also obvious that the effect of axonal activity onto the extracellular medium is additive. In the following, we use the equivalent representation in terms of conductivities and fibre density. The fibre density ρ can be defined as ρ = g −2 A ax /(g −2 A ax + A ex ), where g is the g-ratio of the fibres, i.e. axonal diameter divided by fibre diameter (axon plus myelin). The perturbation of the extracellular medium by a spike in the n th axon is thus given by We note here that we only consider axons with circular shape, therefore A ax can be replaced by the square of the fibre diameter d n .

Perturbation of the membrane potential of a passive axon by a spike in a contiguous one
A spike can be regarded as a travelling wave along the axon. If the axon is sufficiently homogeneous, then this wave appears stationary in the co-moving frame. For simplicity, we assume here that the axons under consideration are either homogeneous, which would be the case for unmyelinated axons, or homogenised in the case of periodically myelinated axons. The term homogenised refers to the technique proposed by Basser (1993), which yields compound variables for the cable parameters that depend on the properties of myelinated and unmyelinated segments. We transform the cable equation, Eq.
(1), into the comoving frame by setting ξ = x − ct, where c is the propagation velocity of the spike, which results in Here, we have made use again of φ i = V +φ e , as well as τ = C m /g m and λ 2 = g ax /g m , and · indicates differentiation with respect to ξ . More specifically, τ is the homogenised time constant, and λ is the homogenised length constant (Basser 1993): where indices denote either myelinated segments with length L, or nodes of Ranvier with length l. The node and myelinspecific parameters were chosen to be τ myel = 0.47 ms, τ node = 0.03 ms, λ myel = 1930 √ ln(1/g)d, and λ node = 55 √ d, with d being the axon diameter. These parameters are the same as in Schmidt and Knösche (2019), and all length scales are measured in μm. For simplicity, we set l/L = 0.01. Since the g-ratio is also held fixed at g = 0.6, the only free parameter is the axon diameter d.
If the perturbation is sufficiently small, it will fail to elicit a spike in the passive axon, and the nonlinear term I HH (V ) can be neglected for simplicity. Thus we arrive at a linear ODE that describes the spatial profile of the perturbation in the passive axon: Here, V p indicates the perturbation of the membrane potential caused by φ e . This is an inhomogeneous differential equation of second order. The homogeneous solution to (9) is found to be where h 1 (ξ ) = exp(ξ/ν + ) and h 2 (ξ ) = exp(−ξ/ν − ) are the fundamental solutions to (9), with Since the perturbations are localised, we require that V p → 0 as ξ → ±∞, and therefore V p,hom = 0. Nevertheless, the fundamental solutions serve to identify the particular solution to (9), which is found by the method of variation of parameters. The particular solution can be posed as which is identical to the full solution to (9), V p (ξ ). Since differentiation will yield only one equation for the two unknowns α 1 (ξ ) and α 2 (ξ ), we may make a further assumption regarding the solution structure: Differentiation of (12) and insertion of the resulting terms into (9) then yields Equations (13) and (14) pose a set of two linear equations that can be solved for α 1 (ξ ) and α 2 (ξ ). Subsequent integration results in and Since the integrands vanish at ±∞, we can pose these indefinite integrals as definite integrals on the interval [ξ, ∞) for α 1 (ξ ), and (−∞, ξ] for α 2 (ξ ) instead. This yields the following integral form for V p (ξ ): An alternative representation is the following convolution integral: with the asymmetric convolution kernel given by The convolution kernel and integral are represented here in absolute space instead of the co-moving frame. In Fig. 1, we visualise the convolution kernel, and compare the perturbation of a passive axon that was obtained numerically, with the convolution of this kernel with the second derivative (curvature) of the spatial spike profile. The relatively small difference between the numerical and analytical result can be explained by nonlinear effects not taken into account in the convolution, and the homogenisation of the axonal parameters.
Thus far, we have considered a single spike that evokes a perturbation in a nearby passive axon [which may be regarded as a 'test axon' as in Goldwyn and Rinzel (2016)]. We now seek to extend the notation to accommodate the perturbations evoked by spikes in multiple axons travelling at different velocities. First, we consider the contribution to the perturbation of axon i by a spike in axon j: where φ e, j is the contribution of the spike in axon j to the extracellular potential, and is the convolution kernel specific to the interaction between axon i and axon j, with Since we consider the linear regime (i.e. I HH (V ) = 0), the perturbations are additive, and the full perturbation and extracellular potential are found to be where A is the set of all active, spike-carrying axons.

Piecewise quadratic approximation of spike profiles
To reduce the computational load, we opt for an approximation of the spike profile that can be adjusted via a shape parameter a 1 . In the time domain, the spike is parameterised as follows: The amplitude is set to V max = 110, and the other parameters can be related to a 1 by imposing smoothness conditions, which results in Thus, the spike has a total duration of T s , and its shape can be varied between a triangular shape at large values of a 1 , and a bell-shape at small values of a 1 . The shape of the spike can then be translated into the spatial domain by calling V (t) = V (ξ/c), which results from the co-moving frame transform.
The piecewise quadratic approximation of the spike profile leads to a piecewise constant curvature, which allows us to compute the perturbation exerted by an active axon on its neighbouring axons analytically. The curvature is given by the second derivative: In the co-moving frame, this translates into The quadratic approximation of the spike profile is shown in Fig. 2, alongside its second derivative and the resulting perturbation profile in a passive axon. The perturbation profile can be computed by inserting Eq. (29) into Eq. (5), and inserting the resulting expression for φ e into Eq. (20). The complete mathematical expression for the perturbation in a passive axon driven by a spike in a second, active one is given by the following expression: where G(ξ ) is a function describing the spatial profile of the perturbation: c Resulting perturbation in a passive axon. Parameters: These expressions can be readily translated into absolute space, which is omitted here for brevity.

Calibration of spike propagation model
The spike propagation model was introduced previously for white-matter fibre bundles (Schmidt et al. 2021). Here, we adapt this model to accommodate the specific properties of peripheral nerve bundles. The major difference lies in how the extracellular potential is generated. While in white-matter fibre bundles one has to consider their radial extent due to their large diameter, peripheral bundles can be regarded as quasi-one-dimensional objects in terms of the parameterisation of the extracellular potential. For instance, if we consider a bundle containing 100 axons with diameters of approximately 1 µm each, the resulting bundle diameter is approximately 10 µm. The electrotonic length constant λ of such axons, however, is approximately 1mm. This means that any effects not covered by the 1D representation of the axon bundle can be neglected, which is in line with the finding of Trayanova et al. (1990), whereby the 1D representation holds if the radius of the fibre bundle (≈ 10 µm in our case, when N = 100) is less than 0.03 times the length of the rising phase of an action potential (>1 mm in our case, cf. Fig. 1a). This implies that even at N = 1000 this approximation would be valid. Furthermore, this holds if the diameters of the axons increase, since both the bundle diameter and the length of rising phase increase linearly with the axon diameter (the latter is due to the linear relationship between axon diameter and spike velocity in myelinated axons). The core concept of the SPM is that a spike can be represented by its position on the axon, and its velocity. We remark here that the number of spikes in the SPM is preserved, i.e. that spikes are neither created ('ectopic spikes') nor annihilated by the ephaptic coupling. The velocity is considered constant along the axon in the unperturbed case. Perturbations of the membrane potential caused by spikes in other, contiguous axons, however, have an effect the propagation velocity, which is modelled in linearised form by: This is based on the assumption that the propagation velocity is determined by a putative spike threshold θ via v = f (θ ), with v 0 = f (V thr ) (Fig. 3a). Eq. (33) is then obtained by linearising f (θ ) around V thr , where γ = − f (V thr ). The perturbation of the membrane potential, V p (x), is computed as shown in the previous section, which requires knowledge of the spike's position. Finally, the position of the spike on the i th axon (each active axon carries one spike) is determined bẏ This means that the perturbation at the location of the spike determines the instantaneous propagation velocity. In general, Eq. (34) has to be solved numerically. For computational reasons, we consider the difference between spike positions to compute v(x): Since the spike position x i refers to the point where the membrane potential first deviates from zero, x thr = √ V thr /a 1 v(x i ) is introduced to mark the position where the membrane potential first reaches V thr . Consequently, the parameters γ and V thr cannot be lumped together and have to be treated separately in the calibration.
In total, the SPM has three free parameters, which are the shape parameter a 1 , and the spike threshold V thr and the coupling parameter γ , which set how strongly the perturbation affects the propagation velocity. The spike threshold also defines the spike position along the spike profile, where the spike is sensitive to perturbations. (We note here that although there are two threshold crossings, we only consider the one along the rising phase, which determines the propagation velocity. The falling phase ('tail') of the spike is due to repolarisation processes, which do not affect the spike velocity.) We first fit the shape parameter a 1 by minimising a cost function, which is defined as the sum of squares of the difference between the model spike profile and the spike profile generated with the biophysical model, see Fig. 3b. Only the rising phase of the spike profile is considered here, since the falling phase does not contribute significantly to the perturbation of the extracellular medium (cf. Fig. 1b). There is a clear minimum of the cost function at a 1 = 740, which yields a near perfect match of the spike profiles along the rising phase (Fig. 3b).
To identify realistic values for γ and V thr , we first generate data for the propagation velocities in N = 10 coupled axons using a biophysically realistic model, which is explained in detail in the next section. We use different parameter values of the fibre density and the fibre diameters. Specifically, the fibre density is chosen at ρ ∈ {0.5 0.6 0.7 0.8 0.9}, and the fibre diameters are chosen to be 1 µm for the smallest axon, and (1 + d)µm for the largest axon, with d ∈ {0.1 0.2 0.3}. Thus, we fit the SPM to N p = 15 parameter combinations p n of ρ and d.
To identify the best-fitting set of parameters, we define the following cost function that is to be minimised: where τ SPM,i and τ BPM,i are the axonal delays of the spike in the i th axon generated with the SPM and the biophys-ical model, respectively. This cost function characterises the model discrepancy for a given set of free parameters P = (γ , V thr ). In Fig. 3, we illustrate the cost function, which shows a clear minimum at γ = 2.785 and V thr = 7.05 mV. The discontinuities in the cost function are due to transitions from asynchrony to synchrony. This result indicates that the perturbation of the membrane potential, V p , affects the propagation velocity quite strongly. We note here that in a previous work on white-matter fibre bundles (Schmidt et al. 2021), the optimal parameters were found to be γ = 6 and V thr = 30 mV. This discrepancy can be explained by the fact that there the extracellular potential φ e was used to compute the perturbation, whereas here we use the perturbed membrane potential. The latter has a smaller amplitude due to axial currents, which is implicitly modelled by the convolution kernel.
To further illustrate the match between the SPM and the biophysical model, we plot spike trajectories generated with the SPM (Fig. 3d) for selected parameter combinations of ρ and d, and compare them with the corresponding trajectories generated with the biophysical model (Fig. 3e).

Biophysical model of spike propagation
The biophysical model is based on Eq. (1), as it provides a detailed description of the nonlinear voltage-gated currents. Specifically, we rewrite Eq. (1) as where V n is the membrane potential of the nth axon, and we insert Eq. (5) for φ e . For numerical purposes, space is discretised into 100 µm long segments. In this description, we distinguish between nodal and internodal segments. Internodal segments do not contain nodes of Ranvier, therefore we set λ = λ myel , τ = τ myel , and I HH (V ) = 0. For nodal segments, we compute λ and τ according to Eq. (8), except that we replace the internodal length L by the length of the segment. Voltage-gated ion channels are only located at the nodal segments, and they obey the following equations: with the dynamics of the gating variables given by Hodgkin and Huxley (1952) The reversal potentials in Eq. (38) are chosen to be e Na = 115 mV, and e K = −12 mV (Brill et al. 1977). While the maximum conductivities were set to g Na = 1200 ms/cm 2 and g K = 90 ms/cm 2 by Brill et al. (1977), we chose g Na = 4800 ms/cm 2 and g K = 720 ms/cm 2 to facilitate faster spike propagation. g m in Eq. (37) was set to 33 ms/cm 2 , which corresponds to τ node = 0.03 ms. The equations were solved using the forward Euler method with t = 0.05μs. To avoid boundary effects, the axons were padded with 200 nodal and internodal segments (20 mm in total) on either side. Neumann boundary conditions were applied to both ends, i.e. V (a) = V (b) = 0, with a = −20 mm and b = 120 mm.

Ephaptic coupling in nerve bundles
In this section, we utilise the SPM to study the effect of fibre heterogeneity on the velocity of, and synchronisation between spikes. A special focus here is on fibre bundles with heterogeneous fibre diameters, which are distributed according to either a shifted, uniform distribution or a shifted alpha distribution. Of particular interest is the interplay between this heterogeneity and the strength of ephaptic coupling as expressed by the fibre density. For simplicity, we focus here on the case of spike volleys that engage all axons in the fibre bundle, and that are completely synchronous initially (i.e. emission times are identical). The novelty here lies in the fact that the SPM reduces the computational effort as compared to the biophysical representation of the axon bundle, which allows us to model the interaction within spike volleys in large bundles. The numerical efficiency allows us to perform parameter screenings to describe the behaviour of the SPM in detail. Specifically, we vary the relative amount of extracellular space, which scales inversely with the amount of ephaptic coupling, and we vary the width of the distribution of axon diameters. The fibre bundles are 100 mm long and we record the axonal delays, which are measured as the time difference between spike initiation and spike termination. We compare two types of distributions: a uniform distribution, and a shifted alpha distribution.

Uniform distribution of fibre diameters
The uniform distribution is set up as follows: where d is the width of the uniform distribution, and (·) is the Heaviside step function. The minimum diameter is fixed at 1 µm, which results in a maximum diameter of (1+ d)µm in this distribution. The mean and standard deviation of this distribution is 1 + d/2 and d/ √ 12, respectively. The coefficient of variation (standard deviation divided by mean) is thus d/ √ 3(2 + d). In Fig. 4a, b, we show the effect of fibre density and diameter distribution on the mean and standard deviation of the axonal delays. At small values of d and large values of ρ, spikes synchronise completely, as indicated by a standard deviation of delays close to zero. This is concurrent with an increase in the mean delay as the fibre density increases. The same synchronisation can be observed at larger values of d, yet a higher fibre density ρ is required to achieve full synchronisation. Before full synchronisation sets in, the standard deviation of delays increases with ρ, which indicates that there is no (complete) synchronisation of the spikes within the volley. Rather, slow spikes form a synchronous cluster which slows down, while fast spikes remain asynchronous. We show a delay density plot in Fig. 4c, which illustrates the delay distribution across ρ for d = 0.1 µm. Here, complete synchronisation occurs at ρ ≈ 0.86, yet already at lower values one can observe clustering and slowing down of spikes. Figure 4d shows 2D histograms of delay distributions and axon diameters at the transition from asynchrony to synchrony. While there is a slight increase in delays with increasing ρ in the asynchronous regime, there is a rapid increase in the delays when full synchronisation sets in.

Fibre diameters from a shifted alpha distribution
To test whether these results depend on the type of diameter distribution, we perform the same analysis for a shifted alpha distribution: The minimum value in this distribution is d = 1, but the maximum diameter depends on the number of axons N . Once more, we find a parameter regime at small values of α and large values of ρ where the spikes completely synchronise (Fig. 5a, b). At larger values of α we find again that synchronisation is only partial, in line with the results obtained for the uniform distribution of diameters. This indicates that the specific type of distribution is not important, rather the width of the distribution is essential for the types of synchronisation (complete or partial) between the spikes. To illustrate the route to synchronisation, we detail the results for α = 0.01 in Fig. 5c. Once more, the spike density plot shows that across a wide range of ρ, synchronisation is partial, and only at around ρ ≈ 0.8 it becomes complete. Again, a slow cluster forms at smaller values of ρ, which is seen in Fig. 5d.

Comparison with biophysical model
We have identified parameters at which the spike volleys undergo complete synchronisation. While we have used the biophysical model to calibrate the SPM, that was done for a relatively small number of axons (N = 10). We now focus on the transition to synchrony with uniformly distributed fibre diameters, and compare the delay distribution between the SPM and the biophysical model for N = 200 axons. Due to the long simulation times for the biophysical model, we focus on three different values of ρ at d = 0.1.
The SPM predicts that at ρ = 0.8 and ρ = 0.85 the spikes do not fully synchronise, yet at ρ = 0.85 one can observe the formation of a slow, synchronised cluster. At ρ = 0.9, the spikes are fully synchronised (Fig. 6a). These results are confirmed by the biophysical model (Fig. 6b), albeit slow clusters in the asynchronous regime seem to be more strongly developed (as indicated by higher spike counts in the bin with highest delays).
It is interesting to note that slow spikes pertaining to small axons synchronise before the entire spike volley synchronises. To explain this effect, we investigate how strongly an axon is perturbed by a spike in a contiguous axon depending on its diameter. Since synchronous spikes experience the hyperpolarising phase of each other's perturbation, we compute the minimum of V p for various diameters in response to a spike in an axon with d = 1 µm. In Fig. 6c, we plot the amplitude of this minimum relative to the minimum in an axon with d = 1 µm. Spikes are generated both with the SPM and the biophysical model, although V p is computed with the analytical approach (cf. Fig. 1b). Our results demonstrate unequivocally that the effect on smaller axons is stronger than on larger axons, which leads to earlier entrainment of slow spikes than of fast spikes.

Discussion
The main contribution of the present study is to illustrate the effect of ephaptic coupling and fibre heterogeneity on the synchronisation of spike volleys in peripheral nerve bundles. The results confirm our initial hypothesis that while ephaptic coupling facilitates synchronisation, fibre heterogeneity counteracts this tendency. This is illustrated by the finding that for increasing fibre heterogeneity, the strength of ephaptic coupling (which increases with fibre density) necessary to produce full synchronisation within a spike volley also increases. This effect is independent of the specific type of fibre distribution chosen. At sufficiently high levels of heterogeneity, spike volleys no longer synchronise, even if the nerve bundle is completely filled with nerve fibre (fibre density ρ = 1). Nevertheless, in the absence of full synchronisation a more subtle phenomenon can be observed, which is the clustering and concurrent slowing down of spikes in small axons, while spikes in large axons do not synchronise. A possible explanation for this is that small axons experience a stronger perturbation of their membrane potential than large axons, cf. Fig. 6c. This clustering appears to be more prominent for the shifted alpha distribution, possibly due to the prominent peak at small axon diameters. To obtain these results for large numbers of axons, we have adapted a spike propagation model (SPM) that was previously devised to simulate spike volleys in white-matter fibre bundles (Schmidt et al. 2021). The core assumption of the SPM is that the propagation of a spike is fully characterised by its intrinsic velocity (determined by structural and electrophysiological parameters of the corresponding nerve fibre), and the extracellular potential generated by spikes in nearby nerve fibres. The SPM therefore represents a much simpler model, without having to solve the nonlinear partial differential equation associated with the full biophysical model. Nonetheless, the SPM contains three free parameters that had to be calibrated using the biophysical model. We chose the scenario of ten ephaptically coupled nerve fibres to generate data for delays with the biophysical model, and the discrepancy between the data and the output of the SPM defined a cost function. Minimising this cost function, we were able to find an unambiguous optimal set of parameters for which the SPM matched best the biophysical results. Regardless of the specific choice of parameters, we surmise that synchronisation and clustering of spike volleys is a universal phenomenon in the SPM.
Some limitations apply to the SPM. One main limitation is that spike volleys travel across the fibre bundle without the possibility of generating (ectopic) spikes, or extinguishing them. In the present study, we only consider spike volleys that engage all axons, therefore only the latter would be a possibility. A mechanism of spike generation and extinction could be incorporated into the SPM by defining conditions involving the perturbation of the membrane potential. More precisely, if the membrane of a passive axon is sufficiently depolarised, then a spike is generated there; conversely, if the membrane is sufficiently hyperpolarised, a spike is extinguished. A previous study, based on a spatially extended FitzHugh-Nagumo model, has demonstrated the possibility of emerging ectopic spikes in ephaptically coupled axons (Sheheitli and Jirsa 2020).
Another limitation is that our study is based on a homogenised formulation of axonal morphology, which discounts the precise location and alignment of nodes of Ranvier. While a previous study has shown that the effect of ephaptic coupling was stronger when nodes are aligned, a staggered arrangement of nodes yielded the same qualitative results (Binczak et al. 2001). Node alignment could be taken into account in the SPM by generating data with the biophysical model for different levels of nodal alignment, and fitting the coupling parameter γ to these values. Alternatively, theoretical considerations as in (Binczak et al. 2001) could serve to modulate the ephaptic coupling strength accordingly.
Our study considers relatively narrow axon diameter distributions to illustrate the effect of ephaptic coupling up to complete synchronisation. Peripheral nerves typically have a wide range of axon diameters, typically 1−10 µm (Ikeda and Oka 2012;Eichel et al. 2020). In the context of our results, wide distributions prevent synchronisation and ensure exact temporal coding based on axon diameter, even if axons are densely packed. However, the clustering regime where only slow spikes synchronise could still be applicable, especially in fibre bundles with a selectivity for large or small axons, thus narrowing the distribution within the bundle. Here, the function of thicker and thinner axons could diverge: thinner axons might transmit convergent, synchronous information, whereas thicker axons might transmit time-critical information.
Finally, we would like to contemplate the possibility to confirm our results experimentally. While it is difficult to record spikes directly in peripheral nerves, it is easy to record the EP generated by nervous activity using surface electrodes. The EP of a spike volley can be computed using the line-source approximation (Holt and Koch 1999;McColgan et al. 2017). A typical EP waveform generated by a single spike is triphasic, with an initial hyperpolarisation, intermediate depolarisation, and final hyperpolarisation. In the case of spike volleys, this wave form is convolved by the spatial distribution of spikes in the volley. As a consequence, highly synchronised volleys will show a strong triphasic profile, whereas distributed volleys show weaker profiles. The clustering that we have observed would then likely result in a multiphasic EP profile, with a fast, lowamplitude response resulting from fast, asynchronous spikes, and a slow, high amplitude response resulting from slow, synchronised spikes. Interestingly, experiments show multiphasic EP profiles, which hints at the existence of such a clustering regime (Parker et al. 2018). Nevertheless, it is also possible that the multiphasic nature of these EPs is due to feedback mechanisms, and it requires more detailed modelling work to disentangle the different mechanisms.