Myelination Increases the Spatial Extent of Analog-Digital Modulation of Synaptic Transmission: A Modeling Study

Analog-digital facilitations (ADFs) have been described in local excitatory brain circuits and correspond to a class of phenomena describing how subthreshold variations of the presynaptic membrane potential influence spike-evoked synaptic transmission. In many brain circuits, ADFs rely on the propagation of somatic membrane potential fluctuations to the presynaptic bouton where they modulate ion channels availability, inducing modifications of the presynaptic spike waveform, the spike-evoked Ca2+ entry, and the transmitter release. Therefore, one major requirement for ADFs to occur is the propagation of subthreshold membrane potential variations from the soma to the presynaptic bouton. To date, reported ADFs space constants are relatively short (250–500 μm) which limits their action to proximal synapses. However, ADFs have been studied either in unmyelinated axons or in juvenile animals in which myelination is incomplete. We examined here the potential gain of ADFs spatial extent caused by myelination using a realistic model of L5 pyramidal cell. Myelination of the axon was found to induce a 3-fold increase in the axonal length constant. As a result, the different forms of ADF were found to display a much longer spatial extent (up to 3,000 μm). In addition, while the internodal length displayed a mild effect, the number of myelin wraps ensheathing the internodes was found to play a critical role in the ADFs spatial extents. We conclude that axonal myelination induces an increase in ADFs spatial extent in our model, thus making ADFs plausible in long-distance connections.

The mechanism underlying d-ADF in cortical networks relies on Kv1 channel inactivation. The subthreshold depolarization leads to inactivation of axonal Kv1 channels which provokes an increase in the presynaptic spike duration (spike broadening), an increase in the spike-evoked Ca 2+ entry in the presynaptic bouton and an enhancement of the transmitter release (reviewed in Debanne et al., 2013;Zbili and Debanne, 2019). h-ADF is a much faster process relying on the recovery of axonal Nav channels from inactivation (or deactivation). A short hyperpolarization of the presynaptic cell (15-200 ms) leads to the deactivation of axonal Nav channels, which provokes an increase in the presynaptic spike amplitude, an increase in the spike-evoked Ca 2+ entry and an enhancement of synaptic transmission. Importantly, for d-ADF and h-ADF to occur at a specific synapse, membrane potential variations of the soma (depolarization or hyperpolarization) have to propagate to the presynaptic bouton to impact the local spike shape (duration or amplitude). Therefore, the spatial extent of d-ADF and h-ADF and the number of synapses impacted by these phenomena are mainly determined by the axonal length constant. By increasing the axonal membrane resistance, myelination could increase the axonal length constant (Castelfranco and Hartline, 2015;Alcami and El Hady, 2019) and therefore expand the number of postsynaptic cells impacted by ADFs. The space constant of analog-digital modulation in unmyelinated axons has been shown to vary between 145 and 430 µm depending of the nature of the subthreshold signal used to induce ADF, the number of branch points, and the cell type (Alle and Geiger, 2006;Sasaki et al., 2012;Bialowas et al., 2015;Rama et al., 2015). In L5 pyramidal neurons which display a myelinated main axon and unmyelinated collaterals, the length constant of the main axon has been evaluated from 417 µm to 1,180 µm Kole et al., 2007;Christie and Jahr, 2009;Cohen et al., 2020). We hypothesized that the development of myelin sheaths in L5 pyramidal neurons, which occurs mainly between P10 and P25 in rats (Battefeld et al., 2019), should increase the axonal length constant and therefore increase the spatial extent of ADFs.
Using a computational approach, we show here that myelination increases the axon length constant by a factor 3 leading to ADFs expression more than 2 mm away from the soma. In addition, the number of myelin layers wrapping the internodes were found to have a critical impact on ADF spatial extent.

Model Morphology
A multi-compartment model of a 36 days-old rat L5 pyramidal neuron was simulated with NEURON 7.6 (see Supplementary Figure S1 for model morphology). The neuronal morphology was taken from a reconstructed neuron by Hay et al. (2011) available on Neuromorpho.org (Ascoli et al., 2007; Neuromorpho ID: NMO_07763; Neuron Name: C080418A-1-SR). The dendritic tree of this neuron was fully reconstructed while the axonal tree was partially reconstructed (up to 1 mm from the soma). The neuron is composed of a dendritic tree, a soma, an axon hillock, an axon initial segment (AIS) and an axonal tree. The axonal tree is composed of the main axon presenting a diameter of 1.14 µm and six axonal collaterals with a diameter of 0.23 µm. The collaterals connect to the main axon at 128.5 µm, 129.6 µm, 300.9 µm, 301.5 µm, 810.1 µm, and 993.9 µm from the soma. We kept unchanged the morphology except for the axon. In fact, in order to observe the spatial extent of ADF on distal synapses, we extended the main axon (total length in our model: 20 mm) and added four distant collaterals branching the main axon 1,907.4, 2,922.4, 3,937.4 and 4,952.4 µm from the soma. Presynaptic sites containing presynaptic Ca 2+ channels were placed every 8 µm into the axon collaterals (Romand et al., 2011) leading to a total number of 1,982 presynaptic sites into the model. The axial intracellular resistivity was fixed to 150 cm in all the model compartments. The membrane capacitance was fixed to 1 µF/cm 2 in all the model compartments except for myelin sheaths (see below). All simulations were run with 6.25 µs time steps and the nominal temperature of the simulation was 37 • C. The measurement locations for subthreshold voltage fluctuations and AP waveform were the nodes of Ranvier of the main axon (or the equivalent places in the case of unmyelinated and hybrid models). The measurement locations for spike evoked Ca 2+ entry and synaptic transmission were the presynaptic sites located in the axon collaterals, leading to an evaluation of d-ADF and h-ADF at these presynaptic sites.
In all the simulations, the densities of the different channels in the dendrites, the soma, the axonal hillock, and the AIS were taken from a previously published model of L5 pyramidal neurons (Hu et al., 2009; see Table 1). These densities were unchanged in the different simulations. In contrast, the channel densities in the axon were modified according to the presence of myelin, the length of the internodes and the number of myelin wraps. These modifications of channel density were made in order to preserve the axonal spike waveform at resting membrane potential measured in the middle of the axon (10 mm away from the soma; Supplementary Figure S2). The preservation of the basal axonal spike waveform in the different simulations was mandatory in order to compare the various conditions we studied. In fact, ADFs depend on spike shape modifications and are highly sensitive to the basal spike waveform. The densities of axonal channels in the different simulations are specified in Table 2.

Myelin Modeling
To model the myelin sheath, one must consider that one passive plasma membrane can be considered as a parallel RC circuit. To simulate the myelin sheath, we used the ''extracellular'' mechanism of the Neuron 7.6 software into the internodes of the model. This mechanism adds a layer of RC circuit to the internodes (Supplementary Figure S3). The passive conductance G my and the capacitance C my of this added layer depend on the number of myelin wraps. We assumed that the passive membrane conductance and the membrane capacitance of the myelin plasma membrane are equal to those of the axon, i.e., G ax , and C ax . As one myelin wrap is composed of two myelin membranes, a myelin sheath composed of one myelin wrap presents a conductance G my = G ax /2 and a capacitance C my = C ax /2. This reasoning extended to a myelin sheath composed of n myelin wraps gives G my = G ax /(2 * n) and C my = C ax /(2 * n). To evaluate the number of myelin wraps we used a g-ratio of 0.698 which has been evaluated in L5 pyramidal neurons of adult rats (Cohen et al., 2020), the diameter of the main axon (1.14 µm), the thickness of the periaxonal space (12.3 nm; Cohen et al., 2020) and the thickness of one plasma membrane (7.5 nm). The g-ratio is the ratio of the inner axonal diameter to the total outer diameter. Therefore, g = r ax r ax +t p +t my , where g is the g-ratio, rax is the axonal radius, tp is the periaxonal space thickness and tmy is the myelin sheath thickness. From this equation, we can see that: t my = r ax * 1 g − 1 − t p . With a g of 0.698, r ax of 0.57 µm and t p of 0.0123 µm, we obtained t my = 0.234 µm. Assuming a thickness of 7.5 nm for one plasma membrane, we obtained that the myelin sheath is constituted of 31.2 plasma membranes. As one myelin wrap is constituted of two plasma membranes, we concluded that the myelin sheath is composed of 15.6 myelin wraps. We chose to apply a value of 15 myelin wraps in our main model of myelinated axon which is in the range of the values observed with electron microscopy (Cohen et al., 2020). Consequently, G my = G ax /30 and C my = C ax /30 in our main model of a myelinated axon. The values of myelin conductance and myelin capacitance in the different simulations are listed in Table 2.

Postsynaptic Responses
To obtain the postsynaptic responses, we used Alpha Synapse Point Processes from Neuron 7.6 inserted into postsynaptic cells. The weights of the synapses were calculated using the charge of the spike-evoked Ca 2+ entry in the presynaptic sites with the following power law: where W is the synaptic weight, A is a scaling factor and Q Ca is the charge of the spike-evoked Ca 2+ current (Scott et al., 2008). Therefore, an increase in the Ca 2+ entry produced by an increase in presynaptic spike amplitude or duration led to an increase in the postsynaptic response amplitude.

Current Injection
AP was produced by a 3 nA current during 3 ms. The 15 mV subthreshold depolarization was produced by a 322 pA current during 10 s. The 15 mV hyperpolarization was produced by a −344 pA current for 200 ms. All these currents were injected into the soma.

Space Constants Calculation
The different phenomena we observed (subthreshold depolarization, subthreshold hyperpolarization, depolarizationinduced AP area increase, hyperpolarization-induced AP overshoot increase) propagated decrementally into the main axon. Due to the presence of axon collaterals and myelinated internodes, the propagations did not follow exactly a monoexponential decay, as it would be expected if the axon was a simple cable. However, in previously published experimental studies, the space constants were evaluated using fit with monoexponential decaying functions. Therefore, in order to compare the values obtained into our model and the values found in previous experimental studies, we chose to define the space constant as the distance at which a given phenomenon reaches 37% of its original value. For example, when we depolarized the soma by 15 mV, we defined the space constant of subthreshold depolarization propagation as the axonal point where the depolarization has decreased to the value of 5.55 mV (15 * 0.37 = 5.55). Importantly, in some conditions, the phenomena did not display a monotonic decay into the axon but evolved in a biphasic manner. We chose to not provide any space constant measurement for these cases.

RESULTS
In order to evaluate the impact of axonal myelination on ADFs spatial extent, we used a computational approach. We used a reconstructed morphology of an L5 pyramidal neuron from a young adult rat (Hay et al., 2011, see ''Materials and Methods'' section for details). The main axon is partially reconstructed (up to 1 mm from the soma) and present six collaterals located respectively at 128.5 µm, 129.6 µm, 300.9 µm, 301.5 µm, 810.1 µm, and 993.9 µm from the soma. In order to observe the propagation of ADFs to long-range connections, we extended the main axon (up to 20 mm from the soma) and we added four distant collaterals located 1,907.4, 2,922.4, 3,937.4 and 4,952.4 micrometers from the soma (the distant collaterals were 1 mm length). Supplemental parts diameters had the same values  as the original axon described in Hay et al. (2011): 1.14 µm for the main axon and 0.23 µm for the collaterals. In our model, dendrites, soma and axonal hillock contained Nav1.2 channels, potassium delay-rectifier channels (Kv) and leak channels. AIS contained Nav1.2 channels, Nav1.6 channels, potassium delayrectifier channels (Kv) and leak channels. The axon contained Nav1.6 channels, Kv1 channels and leak channels. To evaluate the effect of spike waveform on synaptic transmission, the axon collaterals also displayed presynaptic sites every 8 µm (Romand et al., 2011) which, contained a weak density of P/Q type calcium channels (1 pS/µm 2 ). Given the axonal tree extension, the model presented 1,982 presynaptic sites distributed all along the axon collaterals. Spike-evoked Ca 2+ entry and corresponding postsynaptic responses were computed at these presynaptic sites (see ''Materials and Methods'' section). In all our simulations, the densities of channels in the dendrites, the soma, the axon hillock and the AIS were kept unchanged and were taken from a previous model of L5 pyramidal neurons (Hu et al., 2009; Table 1). We only modified the reversal potential of leak channels (E leak = −69.5 mV) to obtain a resting membrane potential of −70 mV. In a first step, we simulated a non-myelinated axon. In this model, the densities of Nav1.6 and Kv1 channels are homogenous all along the axonal tree (g Na1.6 = 370 pS/µm 2 and g Kv1 = 26.8 pS/µm 2 ; Table 2: Non-Myelinated axon). The density of leak channels was homogeneous all along the axonal tree and equal to its value in other compartments (0.333 pS/µm 2 ). However, in order to maintain the resting potential at −70 mV, we had to fix E leak at −38.3 mV in the axon.

Both d-ADF and h-ADF Are Expressed in the Model
The d-ADF occurs when the presynaptic neuron is depolarized before spike generation leading to an increase in transmitter release ( Figure 1A). In pyramidal neurons, d-ADF is a slow process needing several seconds of depolarization to produce an increase in synaptic transmission Kole et al., 2007;Sasaki et al., 2012;Bialowas et al., 2015). This slow time-constant is explained by the slow inactivation time-constant of Kv1 channels (1,500 ms; Shu et al., 2007b). In fact, d-ADF is a Kv1-dependent mechanism in pyramidal neurons: the subthreshold depolarization entails the inactivation of axonal Kv1 channels, leading to the broadening of the axonal spike, an increase in spike-evoked Ca 2+ entry at the presynaptic bouton and an increased transmitter release ( Figure 1A). Importantly, the subthreshold depolarization entails also the inactivation of axonal Nav channels leading to a decrease in the axonal spike amplitude . However, the interplay between the spike broadening and the decrease in spike amplitude leads to an overall increase in spike area leading to an enhancement of presynaptic Ca 2+ entry and transmitter release Kole et al., 2007;Bialowas et al., 2015). To verify the presence of d-ADF in the model, we injected current into the soma in order to produce a spike either generated from the resting membrane potential (V m = −70 mV) or after a 10 s subthreshold depolarization of the soma at −55 mV. Then, we measured the subthreshold A long somatic depolarization provokes the inactivation of axonal Kv1 channels leading to a broadening of the presynaptic spike and an increase in the transmitter release. (B) Schematic description of the h-ADF mechanism. A brief somatic hyperpolarization provokes the deinactivation of axonal Nav channels leading to an increase in both the presynaptic spike amplitude and the transmitter release. (C) d-ADF is present in the model. A spike generated from the depolarized Vm (red trace) is broader than the spike generated from the resting membrane potential (black trace). Therefore, the spike generated from a depolarized Vm leads to a larger presynaptic Ca 2+ entry and a bigger EPSP in the postsynaptic cell. Note that the depolarization is induced by somatic current injection and that the traces were recorded in a presynaptic site, 134 µm from the soma. (D) h-ADF is present in the model. A spike generated from the hyperpolarized Vm (purple trace) presents a larger amplitude than the spike generated from the resting membrane potential (black trace). Therefore, the spike generated from a hyperpolarized Vm leads to a larger presynaptic Ca 2+ entry and a bigger EPSP in the postsynaptic cell. Note that the hyperpolarization is induced by somatic current injection and that the traces were recorded in a presynaptic site, 134 µm from the soma. depolarization value, the spike waveform, the spike-evokedcalcium entry and the postsynaptic response amplitude at the first presynaptic site (located 134 µm from the soma). We found that a 15 mV subthreshold depolarization of the soma propagated decrementally into the axon leading to a value of 9.43 mV at the presynaptic site ( Figure 1C). This depolarization provoked the inactivation of axonal Kv1 channels leading to a broadening of the presynaptic spike, and the inactivation of Nav1.6 channels, producing a decrease in the presynaptic spike amplitude. Overall, the spike area was increased by 15.8%, leading to a 13.4% increase in spike-evoked calcium entry and a 36.8% increase in synaptic transmission ( Figure 1C). Therefore, d-ADF amounted to 136.8% of the control amplitude at 134 µm from the soma in our model, which is similar to previously published d-ADF in pyramidal neurons Kole et al., 2007;Bialowas et al., 2015).
The h-ADF is a much faster process that has been described in CA3 and L5 pyramidal neurons . A 200 ms hyperpolarization of the presynaptic neuron is enough to entail a recovery from inactivation of axonal Nav channels (Nav1.6 in pyramidal neurons), leading to an increase in axonal spike amplitude which produces an increase of spike-evoked Ca 2+ entry and transmitter release ( Figure 1B). Due to the slow recovery from inactivation of Kv1 channels, a 200 ms hyperpolarization has no effect on their availability. Therefore, h-ADF is a purely spikeamplitude dependent phenomenon . To verify the presence of h-ADF in the model, we injected current into the soma in order to produce a spike either generated from resting membrane potential (V m = −70 mV) or after a 200 ms subthreshold hyperpolarization of the soma at −85 mV. Then, we measured the hyperpolarization value, the spike waveform, the spike-evoked-calcium entry and the postsynaptic response amplitude at the first presynaptic site (located 134 µm from the soma). We found that a 15 mV hyperpolarization of the soma propagated detrimentally into the axon leading to a value of 9.13 mV at the presynaptic site ( Figure 1D). This hyperpolarization provoked the recovery from inactivation of axonal Nav1.6 channels, leading to a 17.8% increase in the presynaptic spike overshoot (without modification of its duration), an 11.9% increase in the spike-evoked Ca 2+ entry and a 32.7% increase in synaptic transmission ( Figure 1D). Therefore, h-ADF is similar to what has been obtained experimentally in pyramidal neurons .
We concluded that both d-ADF and h-ADF are reproduced by the model.

d-ADF Spatial Extent Is Increased by Myelination
To quantify the spatial extent of d-ADF, we measured the subthreshold depolarization and the AP area all along the main axon. We found that the depolarization propagated with a space constant of 303.4 µm into the main axon, close to previously published values in L5 pyramidal neurons Kole et al., 2007; Figure 2A, black trace). This led to a space constant for depolarization-induced AP area increase of 338.4 µm into the main axon ( Figure 2B, black trace). Then, in order to observe d-ADF, we computed synaptic transmission in the 1982 presynaptic sites located in axon collaterals (Figure 2C, black trace). In this configuration, we found that 31.4% of the presynaptic sites presented a d-ADF of at least 5% ( Figure 2C, black dots). These presynaptic sites were all located in the proximal collaterals and were on average 307 ± 94.5 µm away from the soma. In order to determine the spatial extent of d-ADF in the model, we plotted the value of d-ADF measured at the first presynaptic site of each collateral as a function of the distance from the soma (Figure 2D, black trace). We defined the spatial extent of d-ADF as the point where this curve crossed the value of 105% for d-ADF. We found that the spatial extent of d-ADF was 836 µm (Figure 2D, black trace). Therefore, d-ADF in an unmyelinated axon is a local phenomenon restricted to proximal synapses (i.e., corresponding to an axonal path smaller than 1 mm; Sasaki et al., 2012).
To evaluate the impact of myelination on the d-ADF spatial extent, we simulated myelinated internodes in the main axon while we left the collaterals unmyelinated. In this model, the main axon was myelinated except at the branching points of the collaterals where 1.5 µm-long nodes of Ranvier were simulated. However, if the distance between two branching points was longer than 200 µm, this part of the main axon has been divided into 100 µm-long internodes with nodes of Ranvier in between. Therefore, the myelinated internodes present a length that varies from 99 to 180 µm (mean value of 101.1 µm), which corresponds to previously published values in L5 pyramidal neurons (Arancibia-Cárcamo et al., 2017). To simulate myelinated internodes, we removed the voltage-gated conductance (Nav1.6 and Kv1) from these compartments, and we added a layer of myelin sheath which conductance and capacitance correspond to 15 myelin wraps (see ''Materials and Methods'' section and Supplementary Figure S3 for myelin modeling). In that case, the myelin sheath conductance was equal to the axonal leak channels conductance divided by a factor 30 and the myelin sheath capacitance was equal to the axonal capacitance divided by 30 (see ''Materials and Methods'' section and Table 2: Myelinated axon 1). These values of myelin membrane conductance and capacitance corresponds to what could be expected for a myelin sheath constituted of 15 myelin wraps, a value that we deduced from a g-ratio of 0.698 found in L5 pyramidal neurons (Cohen et al., 2020), the diameter of the main axon in the reconstructed neuron (1.14 µm) and a periaxonal space of 12.3 nm between the internode's plasma membrane and the myelin sheath (Cohen et al., 2020; see ''Materials and Methods'' section for myelin wraps modeling). Finally, the periaxonal space resistivity was set to 53.7 Ωcm (Cohen et al., 2020). The myelinated internodes are separated by nodes of Ranvier. We simulated nodes of Ranvier by increasing voltage-gated sodium and potassium conductance in order to ensure the spike propagation all along the axon ( Table 2: Myelinated axon 1). Importantly, the exact densities of Nav1.6 and Kv into the nodes of Ranvier were determined in order to preserve the spike waveform at resting membrane potential in the middle of the main axon (10 mm from the soma; Supplementary Figure S2). Finally, the density of the leak channels in the nodes of Ranvier was the same as in the non-myelinated axon.
As expected, we found that myelination increased the velocity of the AP conduction (Supplementary Figure S4). More importantly, we found that the subthreshold depolarization of the soma propagated into the myelinated main axon with a space constant of 906.6 µm (Figure 2A, red trace). Therefore, the myelination entailed a 2.99-fold increase in axonal length constant in our model. In consequence, the myelination increased the spatial extent of the depolarization-induced AP area enhancement (Figure 2B, red trace). Importantly, it was not possible to extract a space constant for the depolarizationinduced AP area enhancement because it did not present a monotonic decay in the myelinated axon (see ''Materials and Methods'' section for space constant calculation). In fact, this parameter evolved in a biphasic manner along the axon: it was stable during the first 650 micrometers then it decreased with the distance from the soma (Figure 2B, red trace). This is easily understandable if one takes into account the two opposite effects of the subthreshold depolarization on the AP shape: it increases its duration and reduces its amplitude. The interplay between these two effects creates the biphasic behavior of the AP area increase in our myelinated model. However, despite the impossibility to calculate a space constantly, our model showed that myelination led to a major increase in the axonal portion affected by the depolarization-induced increase in the AP area. Consequently, we found an increase in the presynaptic sites proportion that presents at least 5% d-ADF (36% instead of 31.4%; Figure 2C, red dots). These presynaptic sites were located both on proximal and distal collaterals and were on average 382.7 ± 318.2 µm away from the soma. Importantly, we found that some presynaptic sites presented d-ADF while they were located more than 2,000 µm away from the soma (Figure 2C, red dots). In fact, the spatial extent of d-ADF was found to be 3,124 µm in this configuration (Figure 2D, red trace). Therefore, the myelination of the main axon allowed the propagation of the somatic subthreshold depolarization to distant collaterals and entailed a major increase in the d-ADF spatial extent.

h-ADF Spatial Extent Is Increased by Myelination
To evaluate the spatial extent of h-ADF, we measured the hyperpolarization and AP overshoot all along the main axon. We also measured the synaptic transmission at presynaptic sites located in axon collaterals. In the unmyelinated model, we found that the hyperpolarization propagated with a space constant of 294.2 µm (Figure 3A, black trace). This leads to a space constant for an AP overshoot increase of 360.3 µm (Figure 3B, black trace). 29.6% of the presynaptic sites displayed an h-ADF of at least 5% ( Figure 3C, black dots). These presynaptic sites were all located in the proximal collaterals and were on average 300.3 ± 91.8 µm away from the soma (Figure 3C, black dots). The h-ADF spatial extent was found to be 815 µm (Figure 3D, black trace).
In the myelinated model, the hyperpolarization propagated into the main axon with a space constant of 858.5 µm (Figure 3A, purple trace). Therefore, the space constant of the propagation of subthreshold hyperpolarization is increased by a factor 2.92 compared to an unmyelinated main axon. This leads to an augmentation of the space constant of hyperpolarizationinduced AP overshoot increase ( Figure 3B, purple trace). Consequently, we found an increase in the proportion of presynaptic sites that displayed at least 5% h-ADF (from 29.6 to 34%; Figure 3C, purple dots). These presynaptic sites were located both on proximal and distal collaterals and were on average 365.5 µm ± 281.24 µm away from the soma. Importantly, as for d-ADF, we found that some presynaptic sites still present h-ADF while they were located more than 2,000 µm away from the soma (Figure 3C, purple dots). In fact, we found that the h-ADF spatial extent was 2,894 µm in the myelinated model ( Figure 3D, purple trace). Therefore, similar to its effect on d-ADF, the myelination induced a major increase in the h-ADF spatial extent.

The Enhanced Spatial Extent of ADFs Is Not Due to Hot Spots of Ion Channels at Nodes of Ranvier
One major difficulty of our approach is the fact that Nav1.6 and Kv1 axonal densities were different in the myelinated and the unmyelinated model ( Table 2). Therefore, the increase in ADFs spatial extent in the myelinated axon could arise either from the presence of myelin sheaths at the internodes or from the increase in ion channels densities into nodes of Ranvier. To determine the key parameter controlling the increase in ADF spatial extent, we performed a simulation with a hybrid main axon presenting the same parameters values than the unmyelinated main axon except at specific hot spots located at the same location than nodes of Ranvier of the myelinated model ( Table 2: Hybrid axon). At these hot spots, the parameters values are identical to those of the nodes of Ranvier of the myelinated model. Therefore, this main axon is a hybrid model that presents nodes of Ranvier-like hot spots separated by portions of non-myelinated axon.
In this hybrid model, we found that the subthreshold depolarization propagated into the main axon with a space constant of 301.6 µm, close to the value found in the unmyelinated main axon (303.4). The depolarization-induced enhancement of the AP area propagated into the main axon with a space constant of 327 µm, close to the value found in the unmyelinated main axon (338.4 µm). We found that 31.2% of the presynaptic sites presented a d-ADF of at least 5% and that they were located on average 306.1 ± 93.9 µm from the soma. These values are close to the values we found for the unmyelinated model (31.4% and 307 ± 94.5 µm). Moreover, d-ADF has a spatial extent of 816 µm in the hybrid axon, close to the value of 836 µm found in the unmyelinated axon. Concerning h-ADF, we found that the hyperpolarization propagated into the main axon with a space constant of 293.4 µm which is close to the value found in the unmyelinated main axon (294.2). The hyperpolarization-induced increase in AP overshoot propagated into the main axon with a space constant of 349 µm, close to the value found in the unmyelinated main axon (360.3 µm). We found that 29.6% of the presynaptic sites presented an h-ADF of at least 5% and that they were located on average 299.5 ± 89.3 µm from the soma. These values are close to the values we found for the unmyelinated model (29.6% and 300.3 ± 91.8 µm). Additionally, h-ADF presented a spatial extent of 808 µm in the hybrid axon, close to the value found in the unmyelinated model (814 µm). Therefore, the spatial extents of d-ADF and h-ADF were similar in the non-myelinated model and in the hybrid model. We conclude that the increase in the ADFs spatial extents in the myelinated model was due to the presence of myelinated internodes and was not an artifact due to high Nav1.6 and Kv1 densities in the nodes of Ranvier.

Effect of Myelination Parameters on ADFs Spatial Extent
Several studies showed that the length of the internodes is variable during development and among neuronal types, which may impact the axonal length constant. In L5 pyramidal neurons, they have been shown to vary between 30 and 150 µm from cell to cell (Arancibia-Cárcamo et al., 2017). In order to observe the effects of internodal length modifications on ADFs spatial extents, we divided all the internodes lengths by 2 (average internodal length: 50.59 µm) or 4 (average internodal length: 25.35 µm) while keeping the distance of collaterals from the soma constant ( Figure 4A). For each condition (internodal length divided by 2 or 4), the axonal densities of Nav1.6 and Kv1 channels were modified in order to keep the axonal spike waveform at the resting membrane potential similar to the one in the non-myelinated axon ( Table 2: Myelinated axon 1, 2, 3; Supplementary Figure S2). Moreover, E leak was slightly modified to maintain the resting membrane potential at −70 mV in the axon ( Table 2: Myelinated axon 1, 2, 3). We found a mild decrease in the propagation of subthreshold depolarization with the internodal length reduction (Figure 4A, red trace). Consequently, the internodal length reduction provoked a small decrease in the space constant for depolarization-induced enhancement of the AP area ( Figure 4B). This led to a limited diminution of the d-ADF spatial extent with internodes shortening (Figure 4D). In fact, d-ADF spatial extent was found to be 2,501 µm, 2,830 µm and 3,124 µm for average internodal lengths of respectively 25.35 µm, 50.59 µm and 101.1 µm (Figure 4D). Similarly, we found a decrease in the axonal length constant for hyperpolarization when the internodal length was reduced ( Figure 4A, purple trace). As a result, hyperpolarization-induced AP overshoot enhancement displayed a shorter space constant (Figure 4C) leading to a decrease in h-ADF spatial extent ( Figure 4E). We concluded that internodal length mildly influences ADFs spatial extents in our model.
The number of myelin wraps increases during development (Looney and Elberger, 1986;Battefeld et al., 2019). Moreover, the number of myelin wraps has been found to vary between 5 and 20 in mature L5 pyramidal neurons (Cohen et al., 2020). We, therefore, explored the effect of myelin wraps number on d-and h-ADF spatial extents. For this, we compared the spatial extents of ADFs into axons which internodes are ensheathed by 0, 2, 5, 10, 15 or 20 myelin wraps ( Figure 5A; see ''Materials and Methods'' section for the modeling of myelin wraps number). We had to modify the densities of Nav1.6 and Kv1 channels in the axon in order to keep the axonal spike waveform and the E leak to maintain the resting membrane potential for each myelin wraps number ( Table 2 and Supplementary Figure S2: Non-myelinated axon and Myelinated axon 1, 4, 5, 6, 7). We found that the subthreshold depolarization space constant for an axon which internodes are ensheathed in 0 myelin wraps was 303.4 µm while it was 531.2 µm, 702.2 µm, 834.3 µm, 906.6 µm and 949.8 µm for axons which internodes are ensheathed respectively into 2, 5, 10, 15 and 20 myelin wraps ( Figure 5A, red trace). The depolarization space constant followed the exponential relationship sp = 313.24 + 636.2 * (1 − e −0.189 * N ) where sp is the depolarization space constant and N the number of myelin wraps. Consequently, the spatial extent of both depolarization-induced AP area enhancement and d-ADF greatly increased with the number of myelin wraps ensheathing internodes (Figures 5B,D). In fact, the spatial extent of d-ADF was found to be 836 µm, 1,827 µm, 2,453 µm, 2,840 µm, 3,124 µm and 3,351 µm when the internodes were ensheathed into respectively 0, 2, 5, 10, 15 or 20 myelin wraps. Similarly, we found that the subthreshold hyperpolarization space constant greatly increased with the number of myelin wraps following the relationship sp = 303.02 + 604.3 * (1 − e −0.175 * N ) ( Figure 5A). Therefore, increasing the number of myelin wraps led to an increase in the spatial extent of both hyperpolarization-induced AP overshoot enhancement and h-ADF (Figures 5C,E). We conclude that the number of myelin wraps ensheathing the internodes is a major determinant of both d-ADF and h-ADF spatial extents, suggesting that in the young adult rat ADFs may extend to more than 2 mm away from the soma.

DISCUSSION
In this study, we show using a computational approach that the axonal myelination may expand both the axonal length constant as well as d-and h-ADF spatial extents by around a factor 3. Furthermore, we show that ADFs spatial extents are critically determined by the number of myelin wraps ensheathing the axon and more modestly by the internodal length. Our work, therefore, suggests that myelinated projection paths such as cortico-striatal or cortico-collicular pathways with an axonal distance ranging between 2 and 3 mm may well express dand h-ADF.

Axonal Length Constant
Axonal length constant depends on several geometrical and electrical parameters as the axonal diameter, the presence of myelin, the number of branch-points or the duration of the voltage shift imposed in the soma. A length constant of ∼2 mm has been estimated in myelinated axons of motoneurons (Gogan  , 1983). In hippocampal axons, the length constant is inversely proportional to the number of branch-points (Sasaki et al., 2012). In L5 pyramidal neurons, length constant is 120 µm for a 10 ms depolarization whereas it can reach 1,000 µm for a 200 ms depolarization (Christie and Jahr, 2009), showing that the propagation of voltage along the axon depends on the frequency of the signal imposed in the somatic compartment.
In L5 pyramidal neurons, the value of the main axon length constant has been estimated between 417 and 1,180 µm for a long subthreshold depolarization (0.2-10 s; Kole et al., 2007;Shu et al., 2007a;Christie and Jahr, 2009;Cohen et al., 2020). Based on our modeling, we propose that this large range of reported values is due to the variability of the axonal myelination in the different studies. In fact, the myelination increases the axonal length constant in our model (Figures 2A,B). Importantly, the study from Shu and colleagues was performed on ferrets, a species in which L5 pyramidal neurons myelination starts distally (<350 µm from the soma; Shu et al., 2007a). This could explain the smaller axonal length constant found in this study (417 µm) compared to other studies performed on rats (553-1,180 µm; Kole et al., 2007;Christie and Jahr, 2009;Cohen et al., 2020) in which the axonal myelination starts just after the axonal initial segment (Battefeld et al., 2019). Moreover, we propose that the large range of axonal length constants found in rat L5 pyramidal neurons is due to the variability of myelination between neurons. In fact, the number of myelin wraps can vary between 5 and 20 in mature L5 pyramidal neurons (Cohen et al., 2020). In our model, the axonal length constant increases with the number of myelin wrap leading to a value of 531 µm for 5 myelin wraps and a value of 949.8 µm for 20 myelin wraps ( Figure 5A). However, while we found a 3-fold effect of myelination on axonal length constant in our simulations, this still needs to be experimentally explored by performing recordings of L5 pyramidal neurons treated by demyelinating drugs such as cuprizone (Hamada et al., 2017).
Importantly, axonal length constant also depends on the axonal membrane resistance. One must consider that somatodendritic membrane resistance decreases during the development (Atkinson and Williams, 2009) which leads to a decrease in the dendritic length constant. If this phenomenon also occurs in axons, we may have overestimated the axonal length constant of mature L5 pyramidal neurons. Nevertheless, the value of 0.333 pS.µm −2 for leak channel density used in our model corresponds to a membrane resistance of 30.03 k .cm 2 which is close to the value of 24.6 k .cm 2 recently estimated by Cohen and colleagues in mature L5 pyramidal neurons (Cohen et al., 2020). Moreover, some studies suggest that the axonal membrane resistance may be much lower than dendritic membrane resistance leading to a large axonal length constant (Dover et al., 2016).

ADF Spatial Extent
The axonal length constant and ADFs spatial extents are biophysically related as ADFs tightly depend on the propagation of somatic subthreshold voltage shifts into the axon leading to the modulation of AP parameters in the presynaptic bouton. The space constant of the depolarization-induced AP broadening underlying d-ADF has been estimated to be near 675 µm in L5 pyramidal cell axons (Kole et al., 2007). In our study, we found that the space constant of the depolarization-induced AP broadening is approximately 338.4 µm for the unmyelinated model but was greatly enhanced in the myelinated model ( Figure 2B). In fact, we showed that increasing the length of internodes or the number of myelin wraps, which occurs during the development of L5 pyramidal neurons (Battefeld et al., 2019), leads to an increase in axonal length constant and ADFs spatial extent (Figures 4, 5).

Myelination Parameters
What is the number of myelin wraps in axons? Myelin wraps in peripheral axons linearly depend on the axon diameter and may vary between 10 and 160 (Arbuthnott et al., 1980;Berthold and Carlstedt, 1982). In the CNS, the number of myelin wraps has been estimated to vary between 5 and 20 wraps (Looney and Elberger, 1986;Bakiri et al., 2011;Harris and Attwell, 2012;Snaidero et al., 2014;Arancibia-Cárcamo et al., 2017;Cohen et al., 2020). Our results suggest that the number of myelin wraps has a critical impact on both the axonal length constant and the ADFs spatial extents (Figure 5). In fact, the space constant for subthreshold depolarization and hyperpolarization went from ∼500 µm with 2 myelin wraps to ∼900 µm for 20 myelin wraps in our model ( Figure 5A). As a consequence, the spatial extent for d-and h-ADF was found to increase in proportion (Figures 5D,E). Therefore, we suggest that the variability of the myelin wraps number found in mature L5 pyramidal neurons could entail a variability of ADFs spatial extent in this cell type. During development, the number of myelin lamellae increases from 0 before birth to maximal levels before sexual maturity (Berthold and Carlstedt, 1982). The myelin thickness also depends on electrical activity (Kaller et al., 2017;Suminaite et al., 2019). Therefore, our model suggests that ADFs spatial extent may increase during the development and vary as a function of neuronal electrical activity.
Regarding the internodal length, we found that this parameter, known to critically determine the conduction velocity (Rushton, 1951;Wu et al., 2012), mildly influenced the axonal length constant and ADFs spatial extents in our model (Figure 4).

Extended ADF in Myelinated Axonal Paths in situ?
Due to a large amount of missing data about axonal physiology, our model presents several unconstrained parameters such as the exact axonal channels density and nature at nodes of Ranvier and internodes or the precise axonal membrane resistance. While it seems likely that the electrical isolation provided by the myelination increases the axonal length constant and ADFs spatial extents, experiments are needed to assess if ADFs can occur at distal synaptic connections (>2 mm). In fact, d-ADF and h-ADF have been mostly studied either in unmyelinated axons or in immature brain circuits in which the myelin is weakly present during early development. As a consequence, the presence of ADFs has never been evaluated in fully myelinated axonal paths. Our data suggest that short projection paths as thalamocortical projection may well express d-and h-ADF. Indeed, in this structure myelination starts rapidly during postnatal development, the length of thalamic axons (∼2 mm) is compatible with ADF and pairedrecording from thalamocortical neurons can be obtained in vitro (Hu and Agmon, 2016) and in vivo (Bruno and Sakmann, 2006). Further experimental work will probably help answer this question.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.