Action potential propagation and synchronisation in myelinated axons

With the advent of advanced MRI techniques it has become possible to study axonal white matter non-invasively and in great detail. Measuring the various parameters of the long-range connections of the brain opens up the possibility to build and refine detailed models of large-scale neuronal activity. One particular challenge is to find a mathematical description of action potential propagation that is sufficiently simple, yet still biologically plausible to model signal transmission across entire axonal fibre bundles. We develop a mathematical framework in which we replace the Hodgkin-Huxley dynamics by a spike-diffuse-spike model with passive sub-threshold dynamics and explicit, threshold-activated ion channel currents. This allows us to study in detail the influence of the various model parameters on the action potential velocity and on the entrainment of action potentials between ephaptically coupled fibres without having to recur to numerical simulations. Specifically, we recover known results regarding the influence of axon diameter, node of Ranvier length and internode length on the velocity of action potentials. Additionally, we find that the velocity depends more strongly on the thickness of the myelin sheath than was suggested by previous theoretical studies. We further explain the slowing down and synchronisation of action potentials in ephaptically coupled fibres by their dynamic interaction. In summary, this study presents a solution to incorporate detailed axonal parameters into a whole-brain modelling framework.


Introduction
Neurons communicate via chemical and electrical signals, and an integral part of this communication is the transmission of action potentials along their axons. The velocity of action potentials is crucial for the right timing in information processing and depends on the dynamics of ion channels studding the axon, but also on its geometrical properties. For instance, the velocity increases approximately linearly with the diameter of myelinated axons [1]. Myelin sheaths around axons are an evolutionary trait in most vertebrates and some invertebrates, which developed independently in several taxa [2]. The presence of a myelin sheath increases the velocity of action potentials by enabling saltatory conduction [3]. Long-term, activitydependent changes in the myelination status of axons are related to learning [4]. The functional role of differentiated myelination is to regulate and synchronise signal transmission across different axonal fibres to enable cognitive function, sensory integration and motor skills [5]. White-matter architecture has also been found to affect the peak frequency of the alpha rhythm [6]. Axons and their supporting cells make up the white matter, which has, for a long time, only been accessible to histological studies [7,8]. With the advent of advanced MRI techniques, some of the geometric parameters of axonal fibre bundles have become accessible to non-invasive methods. Techniques have been proposed to determine the orientation of fibre bundles in the white matter [9] as well as to estimate the distribution of axonal diameters [10], the packing density of axons in a fibre bundle [11,12], and the ratio of the diameters of the axon and the myelin sheath (g-ratio) [13].
First quantitative studies were done by Hursh [14] who established the (approximately) linear relationship between action potential velocity and axonal radius in myelinated axons, and Tasaki [3] who first described saltatory conduction in myelinated axons. Seminal work on ion channel dynamics was later done by Hodgkin and Huxley, establishing the voltage-dependence of ion channel currents [15]. The general result of voltage-dependent gating has been confirmed in vertebrates [16], yet a recent result for mammals suggests that the gating dynamics of sodium channels is faster than described by the original Hodgkin-Huxley model, thereby enabling faster generation and transmission of action potentials [17]. In general, parameters determining channel dynamics differ widely across neuron types [18].
Seminal studies into signal propagation in myelinated axons using computational techniques were done by FitzHugh [19] and Goldman and Albus [20]. Goldman and Albus gave the first computational evidence for the linear increase of the conduction velocity with the radius of the axon, provided that the length of myelinated segments also increases linearly with the axonal radius. The linear relationship is supported by experimental evidence [21], although other studies suggest a slightly nonlinear relationship [22]. More recently, computational studies have investigated the role of the myelin sheath and the relationship between models of different complexity with experimental results [23]. One of the key findings here was that only a myelin sheath with finite capacitance and resistance reproduced experimental results for axonal conduction velocity. Other studies investigated the role of the width of the nodes of Ranvier on signal propagation [24,25], or the effect of ephaptic coupling on signal propagation [26][27][28][29][30][31][32][33].
Most computational studies employ numerical schemes, i.e. they discretise the mathematical problem in space and time and use numerical integration methods to investigate the propagation of action potentials. One problem that arises here is that the spatial discretisation must be relatively coarse to ensure numerical stability, which can be remedied to some extent by advanced numerical methods and computational effort [34]. The other problem, however, cannot be remedied that easily: it is the lack of insight into how the model parameters influence the results, since there is a large number of parameters involved. A way to illustrate parameter dependencies in an efficient manner is to use analytical techniques all the while simplifying the model equations and extracting essential features. Studies that use analytical methods are few and far between [35][36][37][38]; yet it is also worth noting that from a mathematical perspective, myelinated axons are similar to spine-studded dendrites, in the sense that active units are coupled by passive leaky cables. An idea that we pick up from the latter is to simplify the ionic currents crossing the membrane [39,40], there at dendritic spines mediated by neurotransmitters, here at nodes of Ranvier mediated by voltage-gated dynamics.
The goal of this article is to use analytical methods to study the influence of parameters controlling action potential generation, and geometric and electrophysiological parameters of the myelinated axon, on the speed of action potentials. The main focus here is on parameters determining the axonal structure. This will be achieved by replacing the Hodgkin-Huxley dynamics with a spike-diffuse-spike model for action potential generation, i.e. ion currents are released at nodes of Ranvier when the membrane potential reaches a certain threshold. These ion channel currents are considered voltage-independent, but we investigate different forms of currents, ranging from instantaneous currents to currents that incorporate time delays. We also investigate ion currents that closely resemble sodium currents measured experimentally. Our aim is to derive closed-form solutions for the membrane potential along an axon, which yields the relationship of action potential velocity with model parameters.
The specific questions we seek to answer here are the following. First, we query how physiological parameters can be incorporated into our mathematical framework, especially parameters that control the dynamics of the ionic currents. We test if parameters from the literature yield physiologically plausible results for the shape and amplitude of action potentials, and test how the ionic currents from multiple nearby nodes of Ranvier contribute to the formation of action potentials. Secondly, we ask how geometric parameters of an axon affect the transmission speed in a single axon, and how sensitive the transmission speed is to changes in these parameters. We seek to reproduce known results from the literature, such as the dependence of the velocity on axon diameter. We also explore other dependencies, such as on the g-ratio, and other microscopic structural parameters resulting from myelination. We compare the results of our spike-diffuse-spike model with the results from a detailed biophysical model recently used to study the effect of node and internode length on action potential velocity [24]. Thirdly, we investigate how ephaptic coupling affects the transmission speed of action potentials, and what the conditions are for action potentials to synchronise. In particular, we examine how restricted extra-axonal space leads to coupling between two identical axons, and how action potentials travelling through the coupled axons interact.

Results
For the mathematical treatment of action potential propagation along myelinated axons, we consider active elements periodically placed on an infinitely long cable. The latter represents the myelinated axon and is appropriately described as leaky cable, whereas the active elements represent the nodes of Ranvier. In mathematical terms, the governing equation is an inhomogeneous cable equation, which describes the membrane potential V(x, t) of a leaky cable in space x (scalar, longitudinal to the cable) and time t in response to input currents: Here, C m and R m are the (radial) capacitance and resistance of a myelinated fibre, and R c is its axial resistance. The term I chan represents the ion channel currents triggered at nodes of Ranvier. The cable equation (Eq (1)) can be reformulated into by multiplying both sides of (1) with R m . The time constant τ and the cable constant λ are parameters determined by the electrophysiological properties of myelin. We choose these parameters in accordance with experimental results and keep them fixed throughout our analysis, see the Methods section for details. The input currents generated by the ion channel dynamics at the nodes of Ranvier is commonly described by a Hodgkin-Huxley framework. However, the Hodgkin-Huxley equations are a challenge to solve analytically, and in order to proceed with our mathematical treatment we opt for a simplified description using threshold-activated currents with standardised current profiles. We analyse different current profiles, ranging from delta-spikes to combinations of exponentials which give a good approximation of the ion currents observed experimentally. We solve the cable equation for these currents analytically which yields the dynamics of the membrane potential describing the resulting depolarisation / hyperpolarisation along the axon. The linearity of the cable equation in V allows us to describe the response to multiple input currents by the superposition of solutions for single currents. A sketch of the framework is shown in Fig 1.

Ion channel dynamics
The classical Hodgkin-Huxley model is described by a set of nonlinear equations which need to be solved numerically. Over the years, it has seen several modifications and improvements such as the one by Frankenhaeuser and Huxley [16], or the incorporation of additional ion currents [41] given the multitude of ion channel types [42,43]. Also, attempts were made to provide better fits by modifying the exponents of the gating variables [44]. In essence, it is difficult to determine what is the 'right' Hodgkin-Huxley model for specific neuron types. For this reason, it seems prudent to go into the opposite direction and to try to simplify the description of the ion channel dynamics.
Two important contributions into this direction are the one by Fitzhugh [45,46] and Nagumo [47], and the one by Morris and Lecar [48]. They provide a framework in which the slow and the fast variables are lumped and thus yield a two-dimensional reduction of the Hodgkin-Huxley model. The ion currents here are still voltage-dependent.
A crucial simplification towards analytically treatable models is the separation of subthreshold dynamics and spike generation in integrate-and-fire models [49,50]. For instance, in the leaky integrate-and-fire model and the quadratic integrate-and-fire model, the time-tospike can be computed analytically, given initial conditions and a threshold value for the membrane potential. The ion currents are then often modelled as delta-spikes since the ion dynamics is fast in comparison to the (dendritic and somatic) membrane dynamics. The spatial extension of the leaky integrate-and-fire model is the spike-diffuse-spike model, in which activity spreads via passive cables.
Here, we consider four forms of channel current models. All of these have in common that the ion current is initiated after the membrane potential has crossed a threshold V thr , and has a predetermined profile. We denote the four scenarios by the letters A, B, C, and D. In scenario A, the ion channel current is released immediately and instantaneously, i.e.
Here, I 0 denotes the overall ion current, t 0 denotes the time when the membrane potential crosses the threshold, and δ(�) is the delta-distribution, or Dirac's delta. In scenario B, the ion current is also released instantaneously, but with a delay Δ: In scenario C, the ionic current is exponential: Here, τ sp is the decay time, and Θ(�) is the Heaviside step function. With scenario D we aim to approximate the ion currents as measured in mammals such as the rabbit [51] and in the rat [52], which can be described by a superposition of exponential currents: A n exp ðÀ ðt À t 0 Þ=t n ÞYðt À t 0 Þ: ð6Þ The currents entering nearby nodes of Ranvier determine the membrane potential at each node, thus forming an action potential. D: The velocity of an action potential is determined by the distance L between two consecutive nodes, and the time difference t sp it takes to reach a given threshold value.
A sketch of all these scenarios is shown in Fig 2, alongside typical depolarisation curves of the membrane potential.

Current influx and separation
According to Kirchhoff's first law, the channel current that flows into the axon, I chan (t) is counter-balanced by currents flowing axially both ways along the axon, I cable (t), and a radial current that flows back out across the membrane of the node, I node , see Fig 3A for a graphical representation. The ratio of currents that pass along the cable and back across the nodal membrane is determined by the respective resistances: with R λ = R m /λ. Throughout the manuscript, the ratio between I cable and I chan is expressed by Based on experimental findings, we assume that the channel density is constant [52], which implies that the total channel current increases linearly with the node length. This is counterbalanced by the fact that the inverse of the resistance of a node, R À 1 node , also increases linearly with its length. At large node lengths, the current that enters the axon saturates, see Fig 3B. We will examine further below how the node length influences the propagation speed. After the membrane potential V reaches the threshold value V thr , the current I is released. A: The instantaneous current is described by a delta-peak at t 0 , when the threshold value is reached. B: The simplest way to accommodate delays or refractoriness is to introduce a refractory period Δ, after which the instantaneous current is released. C: Exponential current with characteristic time scale τ sp . D: A combination of exponential currents describes a realistic current profile. https://doi.org/10.1371/journal.pcbi.1007004.g002

Influence of nearby nodes
During the propagation of an action potential, ion channel currents are released at multiple nearby nodes that affect the shape and amplitude of the action potential. Because of the linear nature of the cable equation, the effect of multiple input currents can be described by linear superposition: where U describes the depolarisation due to the current at a single node with index n. The internode length L, node length l, cable constant λ and cable constant at a node λ n determine the electrotonic distance between nodes. Node indices n are chosen such that the node with n = 0 is centred at x = 0. Nodes with negative n are the ones the action potential has travelled past, and nodes with positive n are the ones the action potential will travel into. Although we consider infinitely long axons, we cut off the sum at n = −N and n = N for computational feasibility, with N = 10 3 . The action potential is not only shaped by the currents from preceding nodes, but also by currents from subsequent nodes that travel back along the axon. Due to the periodic nature of saltatory conduction, the time difference between any two consecutive nodes is assumed to be the same unknown parameter t sp . The effect of distant nodes is dampened by the fact that in addition to passing along myelinated segments, currents from distant sources also pass by unmyelinated nodes, and therefore further lose amplitude. If nodes are relatively short, the current outflux can be regarded as instantaneous across the node as compared to changes in the current, and the total electrotonic distance between two consecutive nodes (measured in units of λ) is then given by L + lλ/λ n , which is already included in Eq (9). Here, λ n denotes the cable constant at a node. Eq (9) describes the temporal evolution of an action potential in a specific location x. In Fig 4 we dissect an action potential using scenario D for the ion channel model, by colour-coding the depolarisation due to individual nodes. It is apparent that the action potential propagation is a collective process with each node regenerating the action potential by a small fraction.

Velocity of action potentials
We now consider the node at x = 0 (n = 0) to reach the firing threshold V thr at t = 0. The relationship between the firing threshold V thr and the time-to-spike t sp is then given by UðnðL þ ll n =lÞ; nt sp Þ; ð10Þ where we have changed the sign of the summation index, i.e. −n ! n. The choice of x = 0 and t = 0 is without loss of generality. Eq (10) is an implicit equation for t sp , which we solve here numerically using Newton's method. The velocity of an action potential is then given by the physical distance between two consecutive nodes, L + l, and t sp : Here we still assume that the activation process at a node is uniform across its entire length. Since a node represents a short section of unmyelinated axon, we estimate the action potential velocity within a node by the action potential velocity in an unmyelinated axon, v n (see Methods section). The resulting velocity then reads We use Eq (12) throughout the manuscript.

Analytical solutions
In mathematical terms, the depolarisation U resulting from the ion channel current at a single site, is a convolution of the current entering the cable with the Green's function of the homogeneous cable equation G(x, t), which describes the propagation of depolarisation along the myelinated segment: Here, x denotes the distance between the site where the current is injected and the site where the membrane potential is recorded. In the following we present the analytical solutions for all the current types.
Scenario A-Fast current. Since the fast current is described by a delta function, the convolution integral turns into the Green's function up to a prefactor: Here, I 0 = 6.6pA/μm 2 is the amplitude of the input current, R λ = R m /λ, and β is the ratio between the current entering the cable and the channel current, as given by Eq (8). I 0 is chosen such that the amplitude of an action potential is approximately 100mV, with all the other parameters chosen as for scenario D with standard parameters, see Methods section.
Inserting Eq (14) into Eq (9), we obtain the spatio-temporal evolution of an action potential for this scenario: with Θ being the Heaviside step function to ensure causality. The threshold condition (Eq 10) then reads Although this is the simplest scenario, it is not obvious how to invert the r.h.s. of Eq (16) to obtain an explicit expression for t sp . In the Methods section we present a linearisation approach, but it is convenient to solve Eq (16) numerically using Newton's method.
Scenario B-Delayed fast current. The membrane dynamics in scenario B is exactly the same as in scenario A, except for an additional offset Δ: R l bI 0 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 4pðt À DÞ The spatio-temporal evolution of an action potential is now given by R l bI 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 ffi ffi ffi ffi ffi 4pðt À nt sp À DÞ and the threshold condition reads ffi ffiffi t p R l bI 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 4pðnt sp À DÞ q exp À n 2 ðL þ ll n =lÞ 2 t Because multiple nodes contribute to the depolarisation, it is possible to find t sp < Δ.

Scenario C-Exponential current.
Here we have to solve the convolution integral of the cable equation with an exponential function, which yields witht ¼ ðt À 1 À t À 1 c Þ À 1 , I representing the imaginary part of the argument, and erf being the error function. In the Methods section we show how to obtain this solution. Eq (20) thus represents solutions for ion currents with instantaneous onset and exponential decay. Hence, the spatio-temporal evolution of an action potential is expressed by ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi and the threshold condition to determine t sp is ffi ffi ffi ffi ffi ffi ffi Scenario D-Combination of exponentials. The linearity of the cable equation allows us to recur to the solution for scenario C to describe the response to currents described by multiple exponentials. Denoting the solution for one exponential input current with time constant τ s by �ðx; t; t s Þ ¼ e À t=t s R l bI 0 2t ffi ffi ffi ffi ffi tt we express the solution to M superimposed exponential currents by We use this formulation to describe both sodium currents and potassium currents with rising and falling phase. The sodium current is expressed as follows: For simplicity, we focus on the case γ = 1, i.e. the biexponential case. Increasing γ would result in increased initial delays, and therefore lower propagation velocities. The parameter γ also affects the normalisation constant C Na,γ , which ensures that the maximum of I chan,Na is I 0,Na . The potassium current is modelled as throughout the manuscript. In the Methods section we describe how to compute the normalisation constants C Na,γ and C K , and how to convert Eqs (25) and (26) into a sum of exponentials. Hence, the spatio-temporal evolution of an action potential is expressed by ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi t À nt sp p þ i ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi t À nt sp t witht ¼ ðt À 1 À t À 1 s Þ À 1 , and C is the problem-specific normalisation constant. The threshold condition to determine t sp is Anticipating results from the next subsection, we found that scenarios A and C yield velocities that are too fast compared with experimental results. Scenario B allows to adjust the propagation speed by tuning the parameter Δ, yet the shape of the action potential is only determined by the parameters from the cable equation, and thus cannot be adjusted to match experimental results. As it is the most realistic and most flexible model for ion channel currents, we decided to select scenario D to study the sensitivity of the propagation speed to structural parameters.

Sensitivity to parameters
Axon diameter. There is a wide consensus that the propagation velocity in myelinated axons is proportional to the axon diameter. This is mostly due to the fact that both the internode length as well as the electrotonic length constant increase with the diameter. One quantity that does not scale linearly with the axonal diameter is the node length, which determines the amount of current that flows into the axon, as well as setting a correction term for the physical and electrotonic distance between two nodes. We find that the latter introduces a slight nonlinearity at small diameters, although at larger diameters the linear relationship is well preserved, see Fig 5A. In Fig 5A we compare the four ion channel scenarios with experimental results obtained by Boyd and Kalu [53]. Scenario A (instantaneous ion channel current) yields velocities that are about one order of magnitude larger than the experimental results. This suggests that the main bottleneck for faster action potential propagation is indeed ion channel dynamics and their associated delays. Introducing a hard delay with scenario B, we find that we can reproduce the experimentally observed range of velocities. With scenarios C and D we introduce temporally distributed ion channel dynamics. The instantaneous onset and exponential decay of scenario C yields velocities that are slightly faster than experimental results.
In scenario D we explore two sets of parameters. The first set of parameters is obtained by using electrophysiological parameters found in the literature. As it is not obvious how to choose the time constants governing the temporal profile of the ion channel currents, we decided to choose them such that the shape of action potentials of our spike-diffuse-spike model match the shape of action potentials of the biophysical model used by Arancibia-Carcamo et al. [24]. The velocities obtained with this set of parameters fall within the range of experimental results. The second set of parameters is obtained by fitting the model parameters to data generated by the same biophysical model (see Methods). The latter yields velocities slightly below the experimental range, but it matches well the results from the biophysical model.
The present framework also enables us to study unmyelinated axons, in which case the current influx must be adapted, in addition to the physical and electrotonic distance between two neighbouring nodes, which is l and l/λ n , respectively. Since λ n is proportional to , the resulting velocity is also to be expected to scale with ffi ffi ffi d p , see Fig 5B. Making the assumption that the membrane conductivity scales linearly with the ion channel density ρ (ρ is measured relative to the ion channel density of a node), the time constant of the unmyelinated axon scales with τ = τ n /ρ, and the cable constant scales with l ¼ l n = ffi ffi ffi r p . We study different ion channel densities, beginning with the same density as in nodes in the myelinated axon, and then reducing the density to 10% and 2% of the original density. We find that reducing the ion channel density also decreases the propagation velocity. For ρ = 1 we find that the propagation velocity is considerably faster than in myelinated axons at small diameters.
Node and internode length. Two geometric parameters that are not readily accessible to non-invasive MRI techniques are the length of the nodes of Ranvier, and the length of internodes. Here we examine the effect of the node and internode length on the speed of action potentials. We assume that the channel density in a node is constant, which is in agreement with experimental results [52]. The channel current that enters the node is proportional to its length, yet the increase of the node length also means that more of this current flows back across the node rather than entering the internodes. Another effect of the node length is the additional drop-off of the amplitude of axonal currents. Node lengths are known to vary between 1μm and 3μm [24]. The length of internodes is known to increase with the fibre diameter [21,22]. This increase can be understood in light of the fact that the cable constant λ is proportional to the fibre diameter, and therefore increasing the internode length ensures that the ratio L/λ remains at a suitable point for signal transmission.
We restrict the analysis to the activation by sodium currents, since potassium currents are slow and only play a minor role in the initial depolarisation to threshold value. The results are shown graphically for scenario D with standard parameters in Fig 6A, and for parameters fitted to the biophysical model by Arancibia-Carcamo et al. [24] in Fig 6B. Changing the threshold value did have a small effect on the maximum velocity, but did not change the relative dependence on the other parameters.
We find that the propagation velocity varies relatively little with changes in the nodal and internodal length. For scenario D with standard parameters, we find that velocities across the investigated range of parameters are above 70% of the maximum, and for the parameters fitted to the biophysical model the sensitivity is even less. Interestingly, we find that decreasing node length and internode length simultaneously, the velocity increases steadily.
In Fig 6C and 6D we show cross-sections of Fig 6B, and compare these results with numerical results from the cortex model used in [24]. There is a good agreement between our model and the biophysical model, with the biggest discrepancies occurring at short node and internode lengths. We assume that these discrepancies arise due to the fact that the biophysical model only uses 50 nodes, whereas we consider N = 1000 nodes to determine the velocity. In the Methods section, we show that reducing the number of nodes significantly alters the results at short node and internode lengths (Fig 13 in Methods section).
Myelin thickness. The relative thickness of the myelin layer is given by the g-ratio, which is defined as the ratio of inner to outer radius. Hence, a smaller g-ratio indicates a relatively thicker layer of myelin around the axon. In humans, the g-ratio is typically 0.6-0.7, although it is also known to correlate with the axon diameter [54]. In our mathematical framework, the gratio affects the electrotonic length constant λ of the internodes, which scales with ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ln ð1=gÞ p . A classical assumption is that the propagation velocity scales in the same manner [1]. Our results suggest (see Fig 7A) that the velocity depends more strongly on the g-ratio. We therefore generalised this relationship to v = κ(ln(1/g)) α , and find (fitting both κ and α) our results best match α = 0.68 (scenario D with fitted parameters). However, the fitted coefficient α also depends on the ratio of internode length and node length, L/l. We find that α increases monotonically with this ratio (see Fig 7B), and approaches zero when L/l approaches zero. The latter represents the case of an unmyelinated axon.
In Fig 8 we present two-parameter plots of the velocity as function of the g-ratio and axon diameter (Fig 8A), and g-ratio and fibre diameter (Fig 8B). If the axon diameter is held constant, the velocity increases monotonically with decreasing g-ratio. However, if the fibre diameter is held constant, then the velocity saturates at around g = 0.5, because decreasing g at constant fibre diameter means decreasing the axon diameter.

Ephaptic coupling and entrainment
We demonstrate here that it is possible to study the effects of ephaptic coupling on action potential propagation within our framework. We choose two axonal fibres as a simple test case, but more complicated scenarios could also be considered using our analytical approach. Ephaptic coupling occurs due to the resistance and finite size of the extra-cellular space. We follow Reutskiy et al. [31] in considering the axonal fibres to be embedded in a finite sized extra-cellular medium (the space between the axons within an axonal fibre bundle). The resulting cable equation for the n th axon reads C m @ðV n À V e Þ @t ¼ 1 R ax;n @ 2 V n @x 2 À ðV n À V e Þ R m þ I chan n ðtÞ; ð29Þ Parameters: fitted parameters (see Table 1 in Methods section). with V e being the potential of the extra-cellular medium. In the Methods section we describe how to obtain solutions to this set of equations.
We explore solutions to Eq (29) in a number of ways, which are graphically represented in Fig 9. We focus on sodium currents as described by scenario D with standard parameters. First, we study how the coupling could lead to entrainment, i.e. synchronisation of action potentials. To this end, we compare the time courses of V 1 (t) and V 2 (t) in a pair of axons, where an action potential is emitted in the first axon at t = 0, and in the second axon at t = Δt. We then compare the t sp in the neighbouring nodes, and find that for any low threshold values V thr the difference between the t sp is less than Δt, meaning the two action potentials are re-synchronising, see Fig 9A. Next, we asked how the coupling affects the speed of two entrained action potentials. Now we set Δt = 0, in which case V 1 (t) = V 2 (t). We compare the depolarisation curves of the simultanously active axons with when only one axon is active, and find that the voltages rise more slowly if two action potentials are present, thus increasing t sp and decreasing the speed of the two action potentials, see Fig 9B. Thirdly, we considered the case when there is an action potential only in one axon, and computed the voltage in the second, passive axon. We find that the neighbouring axon undergoes a brief spell of hyperpolarisation, with a half-width shorter than that of the action potential. This hyperpolarisation explains why synchronous or near-synchronous pairs of action potentials travel at considerably smaller velocities than single action potentials. The hyperpolarisation is followed by weaker depolarisation.

Discussion
We have developed an analytic framework for the investigation of action potential propagation based on simplified ion currents. Instead of modelling the detailed dynamics of the ion channels and its resulting transmembrane currents, we have adopted a simpler notion by which a threshold value defines the critical voltage for the ion current release. Below that threshold value the membrane dynamics is passive, and once the threshold value is reached the ion current is released in a prescribed fashion regardless of the exact time-course of the voltage before or after. We studied four different scenarios, of which the simplest was described by a deltafunction representing immediate and instantaneous current release. The three other scenarios incorporated delays in different ways, from a shift of the delta function to exponential currents and, lastly, combinations thereof. The latter seemed most appropriate considering experimental results.
The simplified description of the ion currents permitted the use of analytical methods to derive an implicit relationship between model parameters and the time the ion current would depolarise a neighbouring node up to threshold value. This involved the solution of the convolution integral of the ion current with the Green's function of the passive cable equation. From the length of nodes and internodes and the time to threshold value between two consecutive nodes (t sp ) resulted the velocity of the action potential.
We only obtained an implicit relationship between the threshold value V thr and the parameter t sp , which needed to be solved for t sp using root-finding procedures. However, in comparison to full numerical simulations, our scheme still confers a computational advantage, as the computation time is about three orders of magnitude faster than in the biophysical model by Arancibia-Carcamo et al. [24]. In the Methods section we have shown that one can achieve a good approximation by linearising the rising phase of the depolarisation curve. We did not explore this linearisation further, but in future work it might serve as a simple return-map scheme for action potential propagation, in which parameter heterogeneities along the axon could be explored.
We used our scheme to study the shape of action potentials, and we found that the ion currents released at multiple nearby nodes contribute to the shape and amplitude of an action potential. This demonstrates that action potential propagation is a collective process, during which individual nodes replenish the current amplitude without being critical to the success or failure of action potential propagation. Specifically, the rising phase of an action potential is mostly determined by input currents released at backward nodes, whereas the falling phase is determined more prominently by forward nodes (cf. Fig 4).
Our scheme allowed us to perform a detailed analysis of the parameter dependence of the propagation velocity. We recovered previous results for the velocity dependence on the axon diameter, which were an approximately linear relationship with the diameter in myelinated axons, and a square root relationship in unmyelinated axons. Although the node and internode length are not accessible to non-invasive imaging methods, we found it pertinent since a previous study [24] looked into this using numerical simulations. Our scheme confirms their results qualitatively and quantitatively, and performing a more detailed screening of the node length and the internode length revealed that for a wide range the propagation velocity is relatively insensitive to parameter variations.
We also studied the effect of the g-ratio on the propagation velocity, which was stronger than previously reported, as we find that the velocity is proportional to (ln(−g)) α with α � 0.7, whereas the classical assumption was α = 0.5 [1]. Furthermore, we found that α depends on the ratio between node length and internode length, which to the best of our knowledge has not been reported before. Intuitively, changing the thickness of the myelin sheath of relatively short internodes has a smaller effect than changing the myelin thickness around long internodes (relative to the node length).
The main results of our spike-diffuse-spike model were compared with the biophysically detailed model recently presented by Arancibia-Carcamo et al. [24]. The latter uses the Hodgkin-Huxley framework and models the myelin sheath in detail, including periaxonal space and individual myelin layers. To enable the comparison between the two models, we fitted parameters of our spike-diffuse-spike model to output of the biophysical model. In spite of the differences in the model setup, we find that the results of the two models agree well.
The framework developed here also allowed us to study the effect of ephaptic coupling between axons on action potential propagation. We found that the coupling leads to the convergence between sufficiently close action potentials, also known as entrainment. It has been hypothesised that the functional role of entrainment is to re-synchronise spikes of source neurons. We also found that ephaptic coupling leads to a decrease in the propagation speed of two synchronous action potentials. Since the likelihood of two or more action potentials to synchronise in a fibre bundle increases with the firing rate, we hypothesise that a potential effect could be that delays between neuronal populations increase with their firing rate, and thereby enable them to actively modulate delays. In addition, we examined the temporal voltage profile in a passive axon coupled to an axon transmitting an action potential, which led to a brief spell of hyperpolarisation in the passive axon, and subsequent depolarisation. This prompts the question whether this may modulate delays in tightly packed axon bundles without necessarily synchronising action potentials. The three phenomena we report here were all observed by Katz and Schmitt [55] in pairs of unmyelinated axons. Our results predict that the same phenomena occur in pairs (or bundles) of myelinated axons.
There are certain limitations to the framework presented here. First of all, we calibrated the ion currents with data found in the literature. This ignores detailed ion channel dynamics, and it is an open problem how to best match ion currents produced by voltage-gated dynamics with the phenomenological ion currents used in this study. Secondly, we assumed that the axon is periodically myelinated, with constant g-ratio and diameter along the entire axon. The periodicity ensured that the velocity of an action potential can be readily inferred from the time lag between two consecutive nodes. In an aperiodic medium, the threshold times need to be determined for each node separately, resulting in a framework that is computationally more involved. Here it might prove suitable to exploit the linearised expressions for the membrane potential to achieve a good trade-off between accuracy and computational effort. Heterogeneities in the g-ratio or the axon diameter would be harder to resolve, as the corresponding cable equation and its Green's function would contain space-dependent parameters. If individual internodes are homogeneous, then one could probably resort to methods used in [36] to deal with (partially) demyelinated internodes. Thirdly, we studied ephaptic coupling between two identical fibres as a test case. Our framework is capable of dealing with axons of different size too, as well as large numbers of axons. In larger axon bundles, however, it might be necessary to compute the ephaptic coupling from the local field potentials, as the lateral distance between axons may no longer allow for the distance-independent coupling we used here. Nevertheless, it would be interesting to extend our framework to realistic axon bundle morphologies, and test if the predictions we make here, i.e. synchronisation of action potentials and concurrent increase in axonal delay, still hold. If yes, then there may also be the possibility that delays are modulated by the firing rates of neuronal populations.

The cable equation
To model action potential propagation along myelinated axons, we consider a hybrid system of active elements coupled by an infinitely long passive cable. The latter represents the myelinated axon and is appropriately described by the cable equation, whereas the active elements represent the nodes of Ranvier whose dynamics are governed by parametrically reduced, phenomenological dynamics.
In general, a myelinated axon can be described by the following cable equation: where V(x, t) is the trans-membrane potential, I chan (V, t) represents the ionic currents due to the opening of ion channels, and x represents the spatial coordinate longitudinal to the cable. C m and R m are the capacitance and resistance of myelinated segments of the cable. Multiplying both sides of (30) with R m yields where τ = C m R m and l ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi R m =R c p are the time constant and cable constant pertaining to the internodes. All model parameters are listed in Table 1. Cable parameters. The capacitance of a cylindrical capacitor (such as a myelin sheath, or the insulating part of a coaxial cable) can be found by considering the following relationship, with g being the g-ratio, i.e. the ratio between axon diameter and fibre diameter. The parameter � denotes the permittivity of the medium. The radial resistance of the cylinder is given by: The parameter ρ describes the resistivity of the cylindrical medium. Experimental values for the capacitance and radial resistance of a myelinated axon are reported in Goldman and Albus [20], with (taking values from [56] and assuming g = 0.8 in the frog) The values for k 1 and k 2 correspond to the following values for permittivity and resistivity: Finally, the axial resistance per unit length along the inner medium of the cylinder is given by where ρ ax = 110Ocm [20] is the resistivity of the inner-axonal medium, and πd 2 /4 its cross-sectional area.
With these constants at hand, we can now define the parameters of Eq (31): l � 9:65 � 10 2 d ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ln g À 1 p ; We treat the axonal diameter d and the g-ratio g as free parameters, and ρ ax , k 1 and k 2 are treated as constants.
Analytical solution. The inhomogenenous cable equation can be written in compact form: with _ V indicating the time derivative of V, and V@ indicating the second spatial derivative of V. Fourier transformation in x yields an ordinary differential equation of the form, where~indicates the Fourier transformed quantity. The homogeneous part of Eq (40) has the solutionṼ The inhomogeneous solution in t can be found by the method of variation of the constant, which yields the following convolution integral in t: The inverse Fourier transform of Eq (42) then yields the following double convolution integral in x and t: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 2l 2 ðt À sÞ Since we assume the nodes of Ranvier to be discrete sites described by delta functions in x, this integral becomes ultimately a convolution integral in time only.
Thus, we can identify the Green's function of the cable equation (Eq (1)) as Gðx; tÞ ¼ 1 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi This is Green's function representing the time evolution of the voltage in a cable due to an instantaneous, normalised input current at distance x at time t = 0. A graphical representation of G(x, t) is given in Fig 10A for various values of x. We note here that the Green's function contains two time scales. The first is the characteristic time scale of the cable, τ, which indicates the voltage decay across the myelin sheeth. The second time constant is x 2 τ/4λ 2 , which is the time it takes exp(−x 2 τ/4λ 2 t) to reach 1/e � 0.37. This time depends on all cable parameters, and if x/λ < 1 it is significantly faster than τ. Hence, if t � τ, the cable equation can be approximated by Gðx; tÞ ¼ 1 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi or, conversely, in the limit t � τ, it can be approximated by Gðx; tÞ ¼ 1 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi See Fig 10B for a comparison.

Nodal properties
Like the myelinated parts of the axon, the Ranvier nodes are characterised by their electrophysiological properties through the membrane resistance and membrane capacitance, denoted by R n and C n , which result in a characteristic length scale λ n and a characteristic time scale τ n . We use the following values for R n [20] and C n [57]: where R n ¼ g À 1 L , i.e. the inverse leak conductance. With τ n = C n R n we obtain a characteristic time of τ n = 33μs. This value is striking, since typical time constants for neurons at dendrites and the soma range from 10ms to 100ms. This can be explained by the higher density of sodium channels at the nodes of Ranvier than at the soma. As reported in [58], there are approximately 1200 channels per μm 2 at nodal segments, and only about 2.6 channels per μm 2 at the soma. Thus, the ratio of ion channel densities between node and soma is nearly 500. We assume here that the conductance scales linearly with the channel density, which is supported by the fact that the membrane resistance is approximately 10kOcm 2 at the soma.
Current influx and separation. The channel current that flows into the axon, I chan (t) is counter-balanced by currents flowing axially both ways along the axon, I cable (t), and a radial current that flows back out across the membrane of the node, I node : The ratio of currents that pass along the cable and back across the nodal membrane is determined by the respective resistances: where R λ is the longitudinal resistance of the axon, defined by R λ = R m /λ. This relationship yields Hence, with the maximum amplitude of the channel current being I 0 , the maximum amplitude of current entering the cable is βI 0 , where we abbreviate

Approximations and analytical solutions
It is, in general, not possible to find closed-form solutions to the Hodgkin-Huxley model due to the nonlinear dependence of the gating variables on the voltage. We therefore focus here on idealisations of the currents generated by the ion channel dynamics, which is described by a function I chan (t).
In mathematical terms, the depolarisation of the neighbouring node is a convolution of the current entering the cable with the solution of the homogeneous cable equation G(x, t), which describes the propagation of depolarisation along the myelinated axon: In the following we present the mathematical treatment for the scenarios introduced in the Results section, and we focus here on an input current at a single site.
Scenario A-Fast current. The (in mathematical terms) simplest scenario is the one in which the ion current is described by the Dirac delta function: Without loss of generality we set the time of the current, t 0 , to zero. The depolarisation along the cable, and specifically at the neighbouring node at distance x is then given by the Green's function of the cable equation itself: If only one current is injected into the cable, the time t sp when the threshold value V thr is reached is given implicitly by Eq 55 yields an implicit relation for t sp and the model parameters. There is no obvious way of solving 55 for t sp explicitly. One can solve it using Newton's method, and test various parameter dependencies by arc-length continuation. However, we explore here the possibility to derive an approximate solution for t sp , and consequently for the axonal propagation speed v, by linearisation of (55).
A suitable pivot for the linearisation is the inflection point on the rising branch, i.e. € V ¼ 0 and _ V > 0. This ensures that the linearisation around this point is accurate up to order Oðt 2 Þ, and error terms are of order Oðt 3 Þ and higher. It also provides an unambiguous pivot for the linearisation. Differentiating (54) twice yields We multiply all terms by t 4 such that the lowest order term in t is of order zero. Since τ is much larger than the rise time of the depolarisation, we disregard terms of order Oðt 3 Þ and higher. The resulting quadratic equation for the inflection point, t i , yields two positive roots, the smaller of which is 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 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi In the limit of x/λ � 1 we can further simplify this expression to give with z ¼ 1=2 À 1= ffi ffi ffi 6 p . The linear equation for the time-to-spike and the firing threshold is then given by The quantities V(t i ) and _ V ðt i Þ can be approximated to be and A comparison of the full nonlinear solution with the linear approximation is shown in Fig  11A. Scenario B-Delayed fast current. Again we consider a fast current, but one which is emitted with a delay Δ after the membrane potential has reached the threshold value. If we denote by t 0 the time of the threshold crossing, then the ionic current is given by However, by simple linear transformation we may also use t 0 to denote the time of the spike. In this case, a spike will be generated after t sp + Δ in the adjacent node, where t sp is the time to the threshold crossing in the same node, given by Eq (55). The speed of a propagating action potential is then given by neglecting finite transmission speeds at nodes. In the limit of t sp ! 0 we obtain the result which implies that action potentials can never travel faster than (L + l)/Δ. However, if multiple neighbours are taken into account, the velocity can be faster than this estimate. For example, in Fig 5A we show results for this scenario with Δ = 30μs. For an axon diameter of d = 1μm (which corresponds to D � 1.67μm with g = 0.6), we obtain a velocity of about 6m/s, whereas (L + l)/Δ is approximately 3.3m/s (with L = 100μm). Scenario C-Exponential current. At this point, we make the assumption that the channel current rises infinitely fast, and drops off exponentially. In mathematical terms, the currents generated by an action potential at a particular node have the following form: where I 0 denotes the amount of current generated by the channel dynamics, and t 0 denotes the time the spike is generated. The Heaviside step function Θ ensures that I chan (t) = 0 for t < t 0 .
Without loss of generality we set t 0 = 0. The propagated depolarisation is now given by the convolution of the exponential function with the Green's function of the cable equation: Here we uset À 1 ¼ À t À 1 þ t À 1 c . We now briefly sketch how to solve this integral. Disregarding prefactors, the integral I to be solved here is of the form A comparison of the linear approximation with the full nonlinear problem is shown in Fig 11B. Scenario D-Combination of exponentials. Scenario C involved a single exponential function to describe the time course of the channel currents. We now explore more complex time-profiles of channel currents, which can be realised by the sum over M exponential time courses with different amplitudes A s and time constants τ s : In particular, we consider current profiles of the form The normalising factor C ensures that the maximum value of I chan (t) is I 0 , which can be determined experimentally. For the sodium current, we use the current density i Na = 50pA/μm 2 , multiplied by the surface area of the node, throughout the manuscript. This current density yields an amplitude of approximately 100mV for action potentials with standard parameters, although it is twice as high as reported in an experimental study [52]. The reason for the experimental values to be lower might be that for the electrophysiological recordings the axons are severed [59], and ion channels are likely to reorganise and redistribute under such conditions.
Eq 83 can be recast in the form The maximum current is reached at and has the amplitude To construct realistic action potentials, we include both sodium and (fast) potassium channels. The sodium gating dynamics of the original Hodgkin Huxley model are governed by a term m 3 h, where m is the activating gating variable, and h is the inactivating gating variable. Schwarz et al. [60] assume that the dynamics of the resulting ion channel currents can be approximated by with C Na,3 being the normalisation constant. Baranauskas and Martina [17] presented data that best fit the Hodgkin-Huxley model with mh, i.e. a linear relationship with the activating gating variable m. In this case, the activation current in our framework reads I chan;Na ¼ I 0;Na C À 1 Na;1 ð1 À exp ðÀ t=t m ÞÞ exp ðÀ t=t h Þ; with C Na,1 being the normalisation constant for γ = 1. The parameters τ m and τ h represent the time constants of the activation and inactivation of the sodium ion channels. The time with γ, τ 1 , and τ 2 as in Eq 83.

Influence of distant nodes
Action potentials are driven by the ionic currents generated at multiple nodes along the axon. Due to the linear nature of the cable equation, the effect of multiple input currents can be described by linear superposition: where U is the r.h.s. of the respective scenario considered, i.e. U(x, t) describes the depolarisation due to the current at a nearby node. To keep with our previous definition, time is defined by setting t = 0 when the neighbouring node reaches threshold. The relationship between the firing threshold V thr and the time-to-spike t sp is therefore given by UðnL; nt sp Þ: ð98Þ The effect of distant nodes is dampened by the fact that in addition to passing along myelinated segments, currents from distant sources also pass by unmyelinated nodes, and thereby further lose amplitude. Because the distance between two points on the cable is given by L/λ in the cable equation, the added distance due to a node with finite length is l/λ n . Therefore, the physical distance between two consecutive nodes is L + l, and their electrotonic distance is L + (λ/λ n )l in units of λ. This leads to the updated equation for the membrane potential, Eq (9) in the Results section.
As we have shown in Fig 4, the formation of an action potential is a collective process that incorporates ion channel currents from multiple nearby nodes. Throughout the manuscript we set N = 10 3 to ensure all currents are incorporated, although for the standard parameters N = 20 would produce very similar results. However, as we show in Fig 13, reducing N can lead to a considerable reduction of the propagation velocity at short internode lengths.
This framework allows us to describe unmyelinated axons as well. Since the internode length is zero in this case, the node length l is now an arbitrary discretisation of the axon. The membrane potential is now described by where the length constant λ in U needs to be replaced by a length constantl that characterises the electrotonic length of the unmyelinated axon. We introduce a parameter ρ that describes the channel density of the unmyelinated axon relative to the channel density of a node of Ranvier. We assume that the conductivity of the axonal membrane scales linearly with the channel density, which implies that the electrotonic length constant of an unmyelinated axon is l ¼ l n = ffi ffi ffi r p , and its time constant ist ¼ t n =r. The velocity of an action potential is now defined as v = l/t sp .
In addition to the correction terms introduced in Eq (9), we also investigate delays that occur at the nodes due to finite transmission speeds. We assume that action potentials travel with velocities v determined by Eq (9) along myelinated segments, and with velocities v n inferred from Eq (99) at nodes. The corrected velocity is then given by Eq (12) in the Results section.
In the limit of small extra-cellular volume and/or highly resistive extra-cellular medium (R À 1 ex ! 0), the coupling parameter is a ¼ 1= P m R À 1 ax;m . We explore this case in the Results section.

Fitting parameters to biophysical model
In order to compare the spike-diffuse-spike model with the biophysical model presented in [24], we generate data points using the biophysical model for the parameters reported therein for the cortex model, and fit our model parameters to these data points. We define a grid of 3 × 3 data points in L − l-space at L = 27μm, L = 82μm and L = 152μm, and l = 0.5μm, l = 1.5μm and l = 3.5μm. On this grid we determine the action potential velocity of the biophysical model, which is treated as data for the fitting procedure. Next, we use the least squares curve fit as implemented in MATLAB to fit the following eight parameters of the spike-diffuse-spike model to the data: λ, τ, λ n , τ n , τ m , τ h , I 0 , and V thr . We use this fitting procedure because there is no direct correspondence between our model and the biophysical model. The latter implements a Hodgkin-Huxley formalism, as well as a detailed model of the myelin sheath that models each membrane individually and includes periaxonal space. We used the code made available on github by the authors of [24].