A mechanism for neurofilament transport acceleration through nodes of Ranvier

INTRODUCTION

Nerve cells extend long cellular processes called axons and dendrites that form electrical connections with other cells throughout the body, thereby establishing the wiring pattern of the nervous system. Communication along these cellular conduits is achieved by the propagation of action potentials, which are waves of membrane depolarization commonly referred to as nerve impulses. Two fundamental mechanisms by which animals can increase the rate of propagation of nerve impulses along axons are to increase axon diameter or to insulate the axons by myelination (Waxman, 1980; Hartline and Colman, 2007). Myelination is a tight spiral wrapping of an axon by a sheetlike extension of a myelinating glial cell. The myelin sheath along a single axon is arranged in contiguous segments called internodes, each formed by a single myelinating glial cell (Schwann cells in the peripheral nervous system, oligodendrocytes in the central nervous system). Each myelinated internode is separated from the next by a short gap of bare axon known as a node of Ranvier, where ion channels that are responsible for initiation and propagation of the nerve impulse (action potential) are concentrated. By clustering the ion channels at nodes and insulating the axon between nodes, myelinated axons are able to propagate nerve impulses in a saltatory manner in which the depolarization at each node spreads rapidly along the myelinated internode to depolarize the next node (Stämpfli, 1954; Salzer, 2003).

For many decades, it has been known that myelinated axons are constricted locally and abruptly at nodes of Ranvier (Hess and Young, 1952; Berthold, 1978) and that the extent of constriction scales with the internodal axon diameter (Rydmark, 1981; Swärd et al., 1995). Using computational modeling, we and others have shown that these constrictions increase the efficiency of saltatory nerve conduction by decreasing nodal capacitance, thereby reducing the internodal caliber required to achieve a given target conduction velocity as much as threefold (Halter and Clark, 1993; Johnson et al., 2015). We also found that there is an optimum theoretical extent of nodal constriction for any given internodal caliber, and that this matches the extent of constriction observed in animals (Johnson et al., 2015). Thus, nodal constrictions appear to be an evolutionary adaptation that confers significant spatial and metabolic efficiency on myelinated axons.

Because of their length, axons are critically dependent on the intracellular transport of organelles and macromolecules for their growth and survival. This movement is called axonal transport (Brown, 2014). In addition to their electrophysiological significance, described above, nodes of Ranvier have important implications for the mechanisms of axonal transport because they represent potential bottlenecks for the movement of axonally transported cargoes. Some of the most abundant cargoes in the axon are neurofilaments, which are long flexible space-filling protein polymers that function to expand axon caliber (Hoffman, 1995). Neurofilaments move bidirectionally along microtubule tracks in an intermittent manner characterized by short bouts of rapid movement interrupted by pauses of varying duration (Wang et al., 2000; Brown, 2014).

Electron-microscopic studies of axons have shown that the number of neurofilaments declines at nodes of Ranvier as much as 10-fold in the largest axons. This local decrease in neurofilament number does not appear to cause the nodal constrictions, because nodal constrictions are also observed in axons lacking neurofilaments (Perrot et al., 2007). However, the constrictions do appear to have important consequences for neurofilaments, because they create potential bottlenecks for the axonal transport of these polymers (Okamura and Tsukita, 1986). To explore how neurofilaments navigate these potential bottlenecks, we used a fluorescence photoactivation pulse-escape technique to analyze the kinetics of neurofilament transport in contiguous nodes and internodes along myelinated axons of adult mouse peripheral nerves ex vivo (Walker et al., 2019). In this approach, axons expressing neurofilament protein tagged with photoactivatable green fluorescent protein (paGFP) are illuminated with violet light to activate fluorescence in a short axonal window, effectively labeling the neurofilaments at that location. Over time, these fluorescent filaments depart from the activated region by the mechanisms of axonal transport with kinetics that is dictated by the stochastic moving and pausing behavior (Trivedi et al., 2007; Li et al., 2014). We found that neurofilaments accelerate locally in nodes and we proposed that this local acceleration permits these polymers to navigate nodal constrictions, analogously to the increase in the current where a river narrows its banks (Walker et al., 2019). However, the mechanism of this local acceleration was unclear.

An interesting observation that may help explain the acceleration of neurofilaments in nodes of Ranvier is that, while the number of neurofilaments declines locally in nodes, the number of microtubules does not (Tsukita and Ishikawa, 1981; Price et al., 1990, 1993; Reles and Friede, 1991; Hsieh et al., 1994). This suggests that most microtubules course through each node from one internode to the next, whereas most neurofilaments do not (Figure 1A). A consequence of this is that the ratio of neurofilaments to microtubules is much lower in nodes than in internodes. In addition, the overall packing density of these polymers is higher in nodes (Price et al., 1990; Reles and Friede, 1991; Hsieh et al., 1994). Collectively, these factors decrease the average distance between nodal neurofilaments and the nearest microtubule (Figure 1, B and C). Thus, one possible mechanism by which neurofilaments could move faster in nodes is simply that they are closer to, and thus have greater access to their microtubule tracks. Here we use computational modeling to test this hypothesis.

MODEL

General description

We consider that neurofilaments are cargoes of axonal transport that move bidirectionally along microtubule tracks. We model neurofilament movement using an extension of a previous cargo-based model (Jung and Brown, 2009; Li et al., 2012), where each neurofilament cycles stochastically between six kinetic states. In the “on-track” states a, a₀, r, r₀, neurofilaments are associated with their microtubule tracks and exhibit movement in anterograde and retrograde directions (states a and r), interrupted by brief pauses in the resting states a₀ and r₀. During the brief pauses, neurofilaments can reverse direction. The dwell times in these pausing states are of the order of seconds to minutes, resulting in a stop-and-go motion of neurofilaments cycling stochastically between the pausing and moving states (Brown et al., 2005). Using computational modeling, we demonstrated that these four states capture the movement of neurofilaments on short time scales on the order of seconds and minutes but not on longer time scales on the order of hours and days (Brown et al., 2005). To account for this discrepancy, we proposed that neurofilaments in the on-track pausing states a₀ and r₀ can switch to corresponding “off-track” prolonged pausing states a_p and r_p, in which we envisioned that neurofilaments were temporarily disengaged from their microtubule tracks, as if parked on the side of the road. The transitions between the on-track and off-track pausing states were dictated by the rate constants γ_on and γ_off. Subsequently, we confirmed the existence of these distinct pausing states experimentally in cultured neurons and peripheral nerve axons ex vivo using the fluorescence photoactivation pulse-escape technique (Trivedi et al., 2007; Li et al., 2014; Monsma et al., 2014; Walker et al., 2019). The resulting six-state model is depicted in Figure 2.

This model has been instrumental in predicting the kinetics of neurofilament transport in vivo. For example, it revealed that a pulse of radiolabeled neurofilaments forms a Gaussian wave that moves and spreads at rates consistent with the published experimental data, demonstrating that the rapid intermittent movement of neurofilaments observed in cultured neurons can explain the population behavior of these polymers in animals (Brown et al., 2005; Jung and Brown, 2009). This model has also allowed us to gain insight into the kinetics of neurofilament transport in the optic nerve (Li et al., 2012) and the local regulation of neurofilament transport by myelinating glia (Monsma et al., 2014). However, a significant shortcoming of the model is that it offers no insight into the mechanistic basis for any differences in kinetic behavior and does not relate neurofilament content to axon caliber. Here, we address this shortcoming by incorporating features of the cytoskeletal organization of the axon into the six-state kinetic model, allowing us to explore the influence of the proximity of neurofilaments to their microtubule tracks on the transport kinetics.

Neurofilament and microtubule organization

Microtubules and neurofilaments are considered to be linear structures arranged in a parallel array coaligned with the long axis of the axon. For simplicity, we assume that the microtubules are stable tracks that extend the entire length of the axonal domain in our simulations. This assumption is supported by serial section electron microscopy of myelinated axons in mouse peripheral nerve, which have demonstrated that axonal microtubules are very long, with average lengths of 370–760 μm (Tsukita & Ishikawa, 1981). We allow the neurofilaments to move forward and backward along these microtubules and to diffuse laterally, that is, in the radial dimension of the axon, when they are not moving.

Each microtubule has 13 protofilaments and therefore is theoretically capable of supporting 13 lanes of traffic. In reality, however, steric considerations are expected to limit the number of neurofilaments that can engage simultaneously at the same location along a microtubule to a maximum of p = 5 lanes (Figure 1D; see below). The microtubules are assumed to be distributed uniformly throughout the axon, meaning that each microtubule has an equal probability of being at any location within the radial dimension of the axon. In this case, the average distance of a microtubule to its nearest neighbor microtubule is given by (Hertz, 1909), where ρ_MT denotes the density of microtubules, that is, the number of microtubules per μm². The microtubule and neurofilament densities depend on the specific type of axon and, in general, on its cross-sectional area and can be obtained from published morphometric studies. For axons of the mouse sciatic nerve, which are modeled in the present study, the internodal microtubule density is about 10/μm² and is almost independent of axon caliber (Reles and Friede, 1991), resulting in an average nearest distance between microtubules of about 158 nm.

The rate of finding a microtubule track

Inspection of electron micrographs of neurofilament-rich axons reveals that the neurofilaments greatly outnumber the microtubules and that consequently some neurofilaments are adjacent to a microtubule and others are not (Friede and Samorajski, 1970; Price et al., 1988; Reles and Friede, 1991). Thus, the probability of moving cannot be the same for all neurofilaments, and neurofilaments that are not next to a microtubule must move laterally to become on track. In our six-state model, the rate at which an off-track neurofilament moves on track is governed by the rate constant γ_on (Figure 2). In our new model, we consider this process to require a diffusive search in the radial dimension of the axon, and thus γ_on becomes an emergent parameter that depends on the average distance between the neurofilaments and their tracks—that is, the relative cross-sectional densities of these cytoskeletal elements. To define the rate of this radial diffusion, we calculate the mean first passage time for an off-track neurofilament to diffuse the average distance to the nearest microtubule, which is given by T = r²/4D (Redner, 2001), where D denotes the diffusion coefficient of a neurofilament in the axonal cytoplasm and r is the average distance to the nearest microtubule (nearest neighbor distance averaged over all angles), given by (see above). This results in the estimate

(1) which could be considered the maximum on-rate, assuming that when a neurofilament meets a microtubule it will always bind to it. While this approximation neglects the geometric effects of microtubule and neurofilament size, it does capture the dependence of the on-rate on the microtubule density. When many neurofilaments are present and each microtubule can be engaged with only a finite number of neurofilaments p, the binding probability for any given neurofilament is reduced. Denoting by N_MT the total number of microtubules at any location, and by N_on the number of neurofilaments engaged with those microtubules (i.e., the number of neurofilaments in the on-track moving (a, r) and on-track pausing (a₀, r₀) states), the number of available microtubules for binding neurofilaments is smaller and given by N_MT - N_on/p, resulting in a reduced on-rate,

(2) where A denotes the cross-sectional area of the axon.

The diffusion constant for radial neurofilament movement in axons is not known and will depend on the details of the cytoskeletal organization. Xue et al. (2015) estimated a radial diffusion coefficient for neurofilaments using the Einstein relation, D = k_BT/f, which connects the diffusion coefficient D with the viscous drag coefficient f, approximated by the viscous drag coefficient of a rigid cylinder of length l and radius a, f = 4π μ/ln(l/a) (Batchelor, 1970) for movement perpendicular to the axis of the cylinder. The dynamic friction coefficient μ depends on the medium the rigid cylinder moves in. Water has a dynamic friction coefficient of 1 cp resulting in a drag coefficient of 1.1 × 10⁻⁸ kg/s and a neurofilament diffusion coefficient of 0.36 μm²/s at room temperature, assuming a neurofilament radius of 20 nm and a length of 5 μm. This diffusion coefficient results in simulated on-rates that are far higher than the observed on-rates, which are on the order of 10⁻⁵-10⁻³/s. Given the likely entanglement of neurofilaments with other structures and other neurofilaments as they move radially in the axon, as well as the additional dissipation of energy by the large numbers of flexible side arms interacting with the intracellular fluid, the above calculations probably underestimate the drag coefficient. For example, a dynamic friction coefficient of 5000 cp, as suggested for a neurofilament gel (Leterrier and Eyer, 1987), results in a drag coefficient of 5.7 × 10⁻⁵ μm²/s and a diffusion coefficient of 7.2 × 10⁻⁵ μm²/s. Because of this uncertainty, we constrain the diffusion coefficient here to produce on-rates on the order of 10⁻⁴/s (Trivedi et al., 2007; Jung and Brown, 2009; Walker et al., 2019) and neurofilament transport velocities on the order of 0.5 mm/d, corresponding to the velocity of neurofilament transport in mouse ventral root and sciatic nerve motor axons in vivo (Xu and Tung, 2001). This yields a diffusion coefficient of D = 2⋅10⁻⁶ μm²/s.

Predicting axon caliber

Our goal is to build a model that, given a certain number of neurofilaments and microtubules, generates an axon with the appropriate cross-sectional area. However, in addition to neurofilaments and microtubules, there are other structures in axons that occupy space (e.g., membranous organelles) and we need to include that space in our calculations. To this end, we devise a strategy that divides the space occupied by the other structures and reallocates it to the microtubules and neurofilaments so that each of them accounts for the space covered by other structures through effective cross-sectional areas. To determine the effective cross-sectional areas A_NF and A_MT of neurofilaments and microtubules, we require that the sum of all effective cross-sectional areas of all microtubules and neurofilaments adds up to the cross-sectional area A of the axon; that is,

(3) where N_NF denotes the number of neurofilaments and N_MT the number of microtubules. With the densities of neurofilaments and microtubules defined by ρ_NF = N_NF/A and ρ_MT = N_MT/A, we find the relation

(4)

A neurofilament is composed of a backbone with a diameter of about 10 nm and side arms about 15 nm in length that are oriented perpendicular to the backbone and generate a lampbrush-like polymer structure with a diameter of about 10 nm + 2 × 15 nm. A microtubule has an actual diameter of about 25 nm. Using geometric considerations, we assume we can fit five neurofilaments around a microtubule (Lai et al., 2018; Figure 1D), giving the fully occupied microtubule track a diameter of about 25 nm + 2 × 40 nm = 105 nm. Thus, a fully occupied microtubule occupies a 6.89-fold greater cross-sectional area than a single neurofilament. We choose the ratio of the effective cross-sectional areas of neurofilaments to microtubules to reflect the ratio of the actual cross-sectional areas of these two structures, A_MT = 6.89A_NF, yielding, in conjunction with Eq. 4, explicit values for the effective cross-sectional areas:

(5)

The sizes of these effective cross-sectional areas will depend on the specific type of axon being modeled and must be determined from morphometric data for the areal polymer densities in that axon type (see Figure 1, B and C). For the present study, we used the data of Reles and Friede (1991), who analyzed the numbers of microtubules and neurofilaments with respect to axonal cross-sectional area in nodes and internodes of adult mouse sciatic nerve. A linear regression analysis of the data in Figure 5 of that study revealed microtubule and neurofilament densities of 10/μm² and 170/μm² in internodes, and 53/μm² and 209/μm² in nodes of Ranvier. These densities give rise to effective cross-sectional areas of and in the internodes, and and in nodes. These are the effective cross-sectional areas we use in the present study (Table 1).

The concept of the effective cross-sectional area allows us to construct a realistic axon where the numbers of neurofilaments and microtubules are associated with the correct axonal caliber. This allows us to study the effects of changing neurofilament influx and microtubule content, as seen in axons subject to radial growth, as well as to study the effects of a local change of microtubule density, such as the node of Ranvier, which is the subject of the present study.

Implementation of the model

We consider a spatial domain consisting of a one-dimensional 800-μm axonal segment containing N_MT microtubules and N_NF neurofilaments. Each microtubule is considered to extend the entire length of the axonal segment, whereas the neurofilaments are shorter. To track the distribution of neurofilaments along the axon in our simulations, we discretize the axonal domain into bins 1 μm in length. The length of each neurofilament is drawn from a distribution with average length 5.5 μm (minimum 1 μm, maximum 43 μm) obtained in a recent experimental study on cultured neurons (Fenn et al., 2018) and discretized to the nearest integral μm. The neurofilaments are not constrained to align with the bins, so neurofilament segments may occupy part of a 1-μm bin; in this case, we consider the bin to be occupied if the filament extends through at least half of the bin and to be empty if the filament extends through less than half of the bin. At each time point in the simulation, we record the number of neurofilaments, N_NF (n) in each bin n, as well as the number of neurofilaments that are on track in each bin, N_on (n).

We consider a mature axon that is no longer growing (Saitua and Alvarez, 1988) and thus can be assumed to be in steady state on the time course of our simulations. At the start of each simulation, we assign the location and the states of the neurofilaments randomly. However, after the model equilibrates for any given set of parameters, our results are not dependent on these initial assignments. In a previous study using the six-state model, we have shown the existence of such stationary solutions (Li et al., 2014). These solutions were associated with a uniform distribution of neurofilament content along the axon. The difference between that model and the expanded one used in this paper is that the on-rate, γ_on, is now an emergent property that depends on the number of neurofilaments, N_on, in the on-track states (see Eq. 8).

Neurofilaments interact by sharing a number of finite lanes of traffic on the microtubules, which introduces a nonlinearity. While nonlinearities could in principle give rise to new solutions of a different nature, such as wavelike solutions or solutions with a nonuniform distribution of axon caliber along the internodes, empirically, we have not observed such solutions in our computations. For example, increasing the abundance of neurofilaments results in a larger caliber and a reduced velocity, but since there is no spatial variation in these parameters, the caliber remains uniform. We hypothesize that the reason that we do not observe new solutions of a different nature in spite of the nonlinearities is that the nonlinearities remain weak. Specifically, the fraction of lanes on the microtubules actually occupied by neurofilaments remains well below its maximum. In other words, low neurofilament traffic density along the microtubules results in only minor restrictions of neurofilament motility that are not sufficient to generate nonstationary behavior, such as a blockage that would result in unbounded axonal swelling.

Axons contain a certain number of microtubules and neurofilaments packed densely in a certain volume. In myelinated axons, the neurofilaments greatly outnumber the microtubules. Thus, we assume that neurofilaments compete for a finite number of microtubule tracks and we allow the constraint that two neurofilaments cannot occupy the microtubule location within any given lane on a microtubule. The effect of this is that the velocity of neurofilament transport along a microtubule track is slowed in proportion to the density of neurofilament traffic on that track. To implement these rules, we consider that on-track neurofilaments have a target velocity v_a,max if moving in the anterograde direction and a target velocity v_r,max if moving in the retrograde direction, but we allow these speeds to be achieved only if there are sufficient unoccupied lanes at that location along the microtubule in the direction of movement.

For example, consider a neurofilament moving anterogradely along the axon from left to right. We denote the position of the 1-μm bin to the right of its leading end by n_r. If all p lanes on each microtubule in bin n_r are occupied (i.e., if N_on(n_r) ≥ pN_MT), then the velocity is set to v_a = 0 and the filament will not move. Otherwise, the actual velocity v_a of this neurofilament in the anterograde direction is assumed to decrease linearly with the number of other on-track filaments at location n_r, and is given by

(6)

The velocity in the retrograde direction is similarly modulated by the number of on-track filaments ahead of the filament; that is, at the bin to the left of its leading end, n_l,

(7)

Consistent with experimental observations of neurofilament transport in cell culture (Wang et al., 2000; Wang and Brown, 2001; Uchida and Brown, 2004) and the computational modeling of neurofilament transport in vivo (Brown et al., 2005; Jung and Brown, 2009; Li et al., 2012), the reversal rates between the two directions of motion (γ_ar and γ_ra) are assumed to be very small (see Table 1).

Nodes of Ranvier are short spatial domains of about 1 μm length, where the axon is not myelinated (Figure 1A). At nodes of Ranvier, the axon is constricted in area and that constriction extends both proximally and distally under the paranodal loops for a few micrometers. For simplicity, we refer to this as the nodal constriction, though technically it includes both the node and flanking paranodes. For the purposes of the current study, we choose 10 μm as the total constricted length, which is within the range encountered for large axons (Swärd et al., 1995). As described above, in these constricted domains, the density of microtubules and neurofilaments is significantly higher. These changed densities result in smaller effective cross-sectional areas of microtubules and neurofilaments. Combining Eqs. 2 and 3, the expression for the rate γ_on(n) at the nth bin along the axon is

(8) where we use the appropriate values for the effective cross-sectional areas A_MT and A_NF for the nodal domain and the internodal domain.

Since each neurofilament in an axonal cross-section in bin n experiences the same rate constants, including the rate γ_on(n), the number of on-track neurofilaments is proportional to the total number of neurofilaments in that cross-section. As a consequence, the on-rate γ_on (see Eq. 8) and the velocities v_a and v_r (Eqs. 6 and 7) are a function of the ratio of the number of neurofilaments to the number of microtubules. This has the consequence that, for a given neurofilament and microtubule density, the average neurofilament velocity depends only on this ratio, whereas quantities such as the flux and the cross-sectional area scale with the number of neurofilaments in the cross-section.

Simulations

The course of the simulations is summarized as follows. We choose the internodal cross-sectional area A_inter of the axon we would like to consider. We use published morphometric data to determine the numbers of neurofilaments N_NF and microtubules N_MT in an internodal axonal cross-section that is associated with that size and type of axon (Reles and Friede, 1991). This determines the internodal ratio of the number of neurofilaments to the number of microtubules, which is critical for our predictions. We then distribute neurofilaments randomly and uniformly along the axon so that the average number of neurofilaments in each cross-section equals the value N_NF.

We implement the kinetic processes governing the moving and pausing neurofilaments in a stochastic manner and track the location and kinetic state of each neurofilament with time. Since a neurofilament transitions between off-track and on-track states with on-rates that depend on the abundance of neurofilaments, its movement is not independent of its neighbors; the motion of one neurofilament can affect the transition rates of other neurofilaments. Thus, at each time step dt = 1 s we sequentially update each neurofilament’s kinetic state and position (including the total and on-track neurofilament numbers N_NF(n) and N_on(n)) in all affected bins before we proceed to the next time step.

When a neurofilament leaves the distal end of the axonal window it is reinserted at the proximal end; that is, we treat the ends of the axonal segment as a periodic boundary. Although this may seem unrealistic at first glance, it simulates an axonal segment in steady state where the rate of neurofilament entry to the axonal window proximally, j_a(0), is matched by the rate of neurofilament departure distally.

The net flux of neurofilaments, that is, the balance of anterograde and retrograde neurofilament fluxes in bin n, is obtained as

(9)

Here we sum the anterograde and retrograde velocities of each neurofilament in the moving states (given in Eqs. 6 and 7) and extending into the nth bin per bin-length, that is, per 1 μm. Note that v_a0(i ) = v_r0(i ) = v_ap(i ) = v_rp(i ) = 0 for the velocities in the paused on-track and off-track states. The anterograde flux j_a(n) is the flux that would have to be injected proximally into the axon if boundary conditions different from periodic were used. We denote by and the mean net and time-averaged anterograde neurofilament flux at bin n.

The ensemble-averaged velocity of all neurofilaments at a certain location (bin n) is obtained as the sum of the velocities of all neurofilaments residing in bin n divided by the number of neurofilaments at that location,

(10) where v(i) denotes the velocity of neurofilament i that extends into bin n.

Since neurofilaments are long polymers that can span both node and internode, we must decide how to determine the location of each filament in order to determine which axonal domain influences its transport kinetics. Most likely, the most important factor is where along the polymers the motor(s) are bound, but since this is not known, we assume that they are randomly distributed. This means that any bias that may be caused by motors residing proximal or distal to the center will average out. Thus, in our model, we assume the motors to bind to the center of each neurofilament.

Pulse-escape experiments

To simulate a fluorescence photoactivation pulse-escape experiment, we mark all segments of the neurofilaments that lie within a short window to simulate the photoactivation of neurofilaments containing paGFP-tagged neurofilament protein. We then track the total length of fluorescent neurofilament polymer remaining in that activated region over time. For filaments that are partially in the activated region and partially outside of it, we mark only the segments that are within the activated region at the start of the simulation. As the photoactivated neurofilaments leave the activation window, the fluorescence intensity within the activation window declines. To compare the kinetics in nodes and internodes, we choose a window at random along one of the flanking axonal internodes and compare it with a window of identical length within the node. We record the fluorescence decay at intervals of 30 s over 40 min of simulation. As in experiments, the fluorescence in the activation window at each time point is normalized to the fluorescence immediately after photoactivation and plotted against time.

Code availability

The simulations of neurofilament transport through model internode and node segments are run using Matlab (version R2018b). The code is available on Github (Ciocanel, 2019).

RESULTS

Neurofilament transport in internodes

We chose to model axons of adult mouse sciatic nerve because there are published morphometric data on neurofilament and microtubule densities and axon caliber in nodes and internodes of these axons (Reles and Friede, 1991). To test our model, we start by simulating an 800-µm-long internode with a cross-sectional area of 10 μm², which corresponds to an axonal diameter of approximately 3.6 μm. Using the morphometric data reported in Table 1 and Figure 5 in Reles and Friede (1991), we found that such an axon would contain in cross-section about 100 microtubules and 1600 neurofilaments, corresponding to a ratio of neurofilaments to microtubules of 16:1. The neurofilament velocities in both anterograde and retrograde directions were modulated by the availability of open lanes along their microtubule tracks as specified in Eqs. 6 and 7, and similarly the on-rate for binding to these tracks was modulated as described in Eq. 8. We initially placed 100 microtubules and 213,500 neurofilaments with the length distribution reported in Fenn et al. (2018) randomly along the 800-µm axonal window so that the number of neurofilaments per cross-section was about 16 times the number of microtubules. The initial kinetic states of the neurofilaments were assigned randomly, but as explained above, the system rapidly equilibrated, so these initial assignments did not influence the results.

Figure 3A illustrates five different realizations of the neurofilament distribution along the axon in contiguous 1-μm bins 1, 2, and 3 d after the start of each simulation, as well as the averages of those five realizations. The neurofilament content fluctuated spatially, with an average SD of approximately 3% about the mean along a single axon, or one stochastic realization. Figure 3B shows histograms of the binned neurofilament content at the corresponding time points, which reflects fluctuations due to the stochastic and asynchronous movement of these cargoes. The magnitude of the fluctuations in neurofilament content about the mean is determined by the number N_NF of neurofilaments in each cross-section and scales inversely with (unpublished data). While the stochastic fluctuations in neurofilament content along the axons were present at each point in time, there was no significant change in the distributions after 1 d. This indicates that the neurofilaments and their kinetics reached a steady state within this time. The kinetic parameters in the model dictate that the neurofilaments spend the majority of their time in off-track states, as previously validated in our experimental and computational studies (Jung and Brown, 2009; Li et al., 2012). Thus, as shown in Figure 3A, for an average number of neurofilaments per cross-section of approximately 1600, the average number of neurofilaments in the on-track states at steady state was 120 (7.5%).

To further explore the equilibration of the dynamics, we ran the internode simulation for 2 h and tracked average kinetic parameters of neurofilament behavior over this time interval. Figure 4 shows the mean velocity and on-rate of the neurofilaments as functions of time (after the simulation was initialized) during the first hours of this sample internode simulation. We obtained the mean velocities in a specified bin n by calculating the average of the instantaneous velocities of all neurofilaments extending over that bin. Since the internodal domain was kinetically homogeneous, we then averaged over all bins in the axonal window to obtain a mean velocity (and similarly on-rate) for each time. It can be seen that both the velocity and the on-rate stabilized within 20 min of the start of the simulation. Averaging over a longer 10-h internode simulation, we obtained a mean net velocity of 0.42 mm/d in the anterograde direction, which is in the range of published reports of 0.12–0.6 mm/d for neurofilament transport in mouse ventral root and sciatic nerve motor axons obtained using radioisotopic pulse labeling (Xu and Tung, 2001; Jung and Brown, 2009). The mean overall neurofilament on-rate was estimated to be 2.52 × 10⁻⁴ s⁻¹, which is similar to estimates of the on-rate from fitting fluorescence photoactivation pulse-escape data in Trivedi et al. (2007) and Monsma et al. (2014). After equilibration for 1 d, the homogeneous time averages of the net neurofilament flux and the anterograde neurofilament influx (Eq. 9) were obtained as and . This meant that for an axon containing 1600 neurofilaments, the net influx was predicted to be approximately 8 neurofilaments per second. The fast equilibration of the mean velocity and on-rate (in less than an hour) reported in Figure 4 shows that the mean values we report were averaged over sufficiently long time intervals. It is also worth noting that these equilibrium values for velocities, fluxes, and rates are independent of the initial conditions. Our model therefore recapitulates the kinetics of neurofilament transport in internodes and can generate predictions about the flux, which have not been measured experimentally.

Neurofilament transport across nodal constrictions

To investigate how nodes of Ranvier affect neurofilament transport, we simulated an axonal region of the same size (800 μm) but denoted a region in the middle of this domain as a nodal constriction, dividing up the domain into two internodes. Nodes of Ranvier measure about 1 μm in length, but the constricted region is longer because it extends into the paranodal regions flanking the nodes. For the present study, we assumed a constriction 10 μm in length (Swärd et al., 1995), which is similar to experimental measurements in large myelinated axons of adult mouse tibial nerve (Walker et al., 2019).

To determine the magnitude of nodal constriction expected for an internodal cross-sectional area of 10 μm² and 100 microtubules, we assumed that all the internodal microtubules ran through the node and then used the morphometric data in Figure 5 of Reles and Friede (1991) to extract the nodal cross-sectional area corresponding to that number of microtubules. This yielded a nodal cross-sectional area of 1.76 μm², which corresponds to a nodal constriction ratio of A_inter/A_node = 10/1.76 = 5.7. Using the standard deviations of the measured cross-sectional areas in the relevant axon size category provided in Table 1 of Reles and Friede (1991) (±1.41 μm² and ±0.47 μm² for the internodal and nodal cross-sectional areas, respectively), the SD of this constriction ratio was estimated to be ±1.3 (Stuart and Ord, 1994). To extract the neurofilament content ratio corresponding to a constriction ratio of 5.7 ± 1.3, we used the morphometric data in Table 1 of Reles and Friede (1991), which yields a mean ratio and SD of 3.8 ± 2.3. This represented the target neurofilament content ratio (ratio of neurofilament number in internode versus node) for our simulations.

As mentioned above, in addition to a decrease in the neurofilament-to-microtubule ratio in nodes, there was also an increase in the neurofilament and microtubule packing densities (Price et al., 1990; Reles and Friede, 1991; Hsieh et al., 1994). This resulted in smaller effective cross-sectional areas for these polymers (see Eq. 5; Table 1). Therefore, the nodal constriction was characterized in our simulations by a higher probability of neurofilaments engaging with microtubules, that is, a higher on-rate γ_on.

The simulation protocol in the presence of a nodal constriction was identical to the protocol we used for the internodal simulations. We initially distributed 213,500 neurofilaments uniformly along the entire 800-μm-long axon with the same length distribution as in the internodal case and assigned the filaments randomly to initial kinetic states. We then adjusted the effective cross-sectional areas of neurofilaments and microtubules in the node as described in Model (see Table 1) to reflect the higher density of these polymers in the nodal domain. The axonal cross-sectional area was governed by the effective cross-sectional areas of these polymers in our model, so that the cross-sectional area of the node decreased slightly. Since we defined the on-rate in terms of the proximity of the neurofilaments to their microtubule tracks, this higher packing density resulted in a higher on-rate (see Eq. 8), and thus neurofilaments left this domain more rapidly than they entered. As a consequence of this imbalance, the nodal neurofilament content declined, resulting in a further decline in axonal cross-sectional area, an even larger microtubule density, and a further increase in the neurofilament velocity. This positive feedback loop resulted in unstable dynamics that continued until the nodal cross-sectional area and neurofilament content had declined enough that the neurofilament fluxes in the internode matched those in the node. Within 1 d, the system attained an equilibrium resulting in a stable nodal constriction. Importantly, this mechanism of generating the node in our simulations was not meant to reflect the physiological mechanism by which they arise in development; it was simply a means to attain the stable constriction that is the focus of this study.

In Figure 5, we illustrate the neurofilament content and cross-sectional area across the nodes at three time points to show the stability and fluctuations in these measures once the system has reached steady state. It can be seen that there was a marked decrease in neurofilament number in the nodal constrictions (Figure 5A), which correlated with a reduction in axonal cross-sectional area (Figure 5B). The predicted ratio of internodal to nodal cross-sectional area (the nodal constriction ratio) was 4.3, which is slightly lower than the value of 5.7 ± 1.3 obtained from the experimental data of Reles and Friede (1991). The number of on-track neurofilaments remained constant across the node (green lines in Figure 5A), reflecting the continuity of the flow of neurofilaments. This means that the decline in neurofilament number in the node was due entirely to a decrease in the number of off-track neurofilaments. In other words neurofilaments move faster across the nodes by spending less time off track. This observation is consistent with findings of pulse-escape experiments in Walker et al. (2019).

It is notable that the decrease in neurofilament number was less abrupt (shallower) than the decrease in cross-sectional area. This is because the axonal cross-sectional area is determined by the effective cross-sectional areas of the polymers, which are modulated abruptly within the nodal domain, whereas the neurofilament content is determined by the neurofilament polymers, which form a staggered overlapping array, with single neurofilaments often spanning the boundary between the nodal and internodal domains. We provide an animation of the nodal simulation at steady state, with neurofilament content recorded at intervals of 30 min, in the Supplemental Video.

Neurofilament transport velocity in nodes and internodes

We further explored neurofilament transport across the nodes of Ranvier in these simulations by calculating average parameters of the dynamics, such as velocities and on-rates. At each time point, we looped through all the 1-μm bins along the entire axonal domain and calculate the average of the velocities and on-rates of the neurofilaments that extend into each bin. This means that for each bin, we considered not only those neurofilaments whose center was occupying that bin, but also neurofilaments in that bin that were centered in other bins. Using this approach, we obtained a mean neurofilament velocity and on-rate in the internodal regions flanking the node of 0.42 mm/d and 2.52 × 10⁻⁴ s⁻¹, respectively, which are identical to the values in the earlier internode-only simulations. In contrast, the mean velocity and on-rate averaged across the nodal constriction were 0.93 mm/d and 6.59 × 10⁻⁴ s⁻¹, respectively, reflecting the improved access of nodal neurofilaments to their microtubule tracks (see Figure 6).

To quantify the effect of nodal constriction on neurofilament content, we calculated the ratio of the mean neurofilament content in the internodes (measured as covering the distances 0–300 and 500–800 μm along the axon to avoid the node and flanking regions) versus the neurofilament content in the center (middle bin) of the node (Figure 7A). We refer to this as the content ratio. Figure 7B shows that this ratio increased over a period of several hours after the start of the simulation and then remained stable (with some fluctuations) around 2.5–3 over a period of 3 d. Note that this is not intended to represent how the neurofilament content ratio develops in vivo, but rather to capture the stability of the node at steady state. The numerically predicted neurofilament content ratio is consistent with the content ratio extracted from Reles and Friede (1991), estimated to be 3.8 ± 2.3 for a 10-μm² axon (see above). We note that this neurofilament content ratio is less than the ratio of internodal to nodal cross-sectional area calculated above (the nodal constriction ratio), because microtubules and neurofilaments are packed more densely in the nodal domain.

For these simulations, we used a neurofilament length distribution obtained from cultured neurons (Fenn et al., 2018), because there are no published measurements of neurofilament length distribution in vivo. Since the neurofilament length distribution may be different in myelinated axons in vivo, we used our model to explore the dependence of the nodal morphometry and kinetics on neurofilament length. As for the data in cultured neurons, we assumed an exponential length distribution, but we varied the average length of this distribution over the range 5–20 μm. We implemented simulations on an 800-μm domain with a 10-μm node as described above. Figure 8, A and B, shows the neurofilament content at day 3 for each length distribution. As expected, the decrease in neurofilament content at the node was less pronounced for neurofilament lengths that exceed the length of the nodal constriction, and more pronounced at shorter neurofilament lengths. Figure 8C shows that there was a corresponding reduction in the neurofilament content ratio with increasing neurofilament length. We also explored the nodal “sharpness” by calculating the ratio D/L in Figure 8D, where D is the depth of the node and L is its length for each simulation (see the cartoon in Figure 7A). Nodal sharpness also decreased with increasing average neurofilament length, making the neurofilament content profiles wider and more shallow. This suggests that the length distribution of neurofilaments in vivo may have significance for the neurofilament distribution across nodes, though we recognize that there are other geometric and cytoarchitectural factors that may also come into play, which are beyond the scope of the current model.

In our simulations, we assumed an axon with a cross-sectional area of 10 μm² and a ratio of neurofilaments to microtubules in the internodes of about 16. This ratio, however, can vary between axons of different calibers and axons of different types of neurons. Table 1 in Reles and Friede (1991) reported that the ratio of neurofilaments to microtubules increases with increasing internodal caliber. For example, for axons with a diameter smaller than 1.5 μm, the ratio of the number of neurofilaments to microtubules in internodes was about 7:1, while for axons with a diameter between 3.6 and 4.2 μm, the ratio was about 16:1. Similarly, Reles and Friede (1991) reported that the nodal constriction ratio also increases with increasing internodal caliber, from a value of about 1.3:1 for axons of diameter smaller than 1.5 μm to a value of 5.6:1 for axons with a diameter between 3.6 and 4.2 μm. This suggests that the nodal constriction ratio increased with increasing ratios of neurofilaments to microtubules. This trend is consistent with our model predictions shown in Figure 9, which demonstrate the impact of increasing the ratio of neurofilaments to microtubules on the neurofilament distribution across the node. Increasing the ratio of neurofilaments to microtubules led to a larger decrease in neurofilament content across the nodal constriction.

Simulated pulse-escape experiments

To test for overall consistency, we used the model to simulate fluorescence photoactivation pulse-escape experiments in contiguous nodes and internodes, mimicking the original experimental approach that we use to demonstrate the acceleration of nodal neurofilaments in vivo (Walker et al., 2019). Specifically, we simulated 10 axons with 5-μm activation windows and track the fluorescence decay from these windows, both inside the node and in an adjacent internode. Figure 10A illustrates the pulse-escape strategy, in which violet light is used to activate the fluorescence of neurofilaments within a short axonal segment (either in the node or in the flanking internode) and then fluorescent filaments depart the activated region over time due to their rapid intermittent movement. In multiple independent studies, we have consistently observed that the fluorescence in the activated region decays with biphasic kinetics, with an initial more rapidly declining phase followed by a transition to a more slowly declining phase (Trivedi et al., 2007; Alami et al., 2009; Monsma et al., 2014; Walker et al., 2019). In our model of neurofilament transport, the initial phase represents the departure of neurofilaments that are on track at the time of photoactivation and thus depart within minutes. After those filaments have cleared the activation window, the decay transitions to a slower phase, which represents the mobilization of filaments that are pausing off track and must move on track before they can depart. Thus, the initial slope of the decay curve is dictated primarily by the ratio of the on-track rate constants, γ₀₁/γ₁₀ (which determines the fraction of neurofilaments in the on-track states in our model), and the slope at later times is dictated largely by the on-rate γ_on (Li et al., 2014). Figure 10B compares the simulated kinetics with the experimental kinetics for myelinated axons of similar internodal caliber (average = 4.4 µm) in adult mouse tibial nerve (Walker et al., 2019). The predicted fluorescence decay was similar to the experimental decay for both nodes and internodes, and consistently within one SD of the experimental averages. As in Walker et al. (2019), the fluorescent neurofilaments in the model left the activated regions faster in nodes than in internodes, as evidenced by the initial slope of decay. Thus, the model can explain both the published morphometric and kinetic data on neurofilament distribution and transport across axonal constrictions at nodes of Ranvier.

DISCUSSION

We have developed a model to test the hypothesis that the local acceleration of neurofilaments in nodes of Ranvier can be explained by a local difference in the access of these polymers to their microtubule tracks. The rationale for this hypothesis was based on two key observations in the published literature: 1) neurofilaments and microtubules are packed more densely in nodes, and 2) the ratio of neurofilaments to microtubules is lower in nodes, due largely to a local decrease in neurofilament number. Our model is based on the simple constraint that a neurofilament must be next to a microtubule in order to move along it and the observation that this is often not the case in neurofilament-rich axons, where neurofilaments greatly outnumber microtubules. On the basis of prior kinetic analyses, we consider that the neurofilaments alternate between distinct kinetic states, which we term on-track and off-track. Neurofilaments in the on-track state move rapidly and intermittently along microtubules, pausing only briefly between bouts of movement until they disengage and become off-track. Neurofilaments in the off-track state exhibit extended pauses while they execute a radial diffusive search for another microtubule, whereupon they engage with that microtubule and move back on track. The average search time for a neurofilament to find a microtubule, which depends on their relative proximity, determines the on-rate. To relate axon caliber to the number of neurofilaments and microtubules, we assigned these polymers effective cross-sectional areas that we extracted from published morphometric data.

To simulate a node, we locally increased the neurofilament and microtubule density in a short segment of the axon by reducing the effective cross-sectional areas of these cytoskeletal polymers according to published measurements. This perturbation resulted in a local increase in the neurofilament on-rate and a local increase in the transport velocity, yielding an emerging decline in the neurofilament content and cross-sectional area of the node, that is, a nodal constriction. The number of filaments that were on track in the node at any point in time was similar to that in the internode (consistent with the fact that the microtubule number does not decrease in the node), but the number of filaments that were off track was reduced, resulting in an increase in the average time spent on track and thus in the average velocity, consistent with the acceleration in neurofilament transport observed experimentally in mouse tibial nerve ex vivo (Walker et al., 2019). The decrease in neurofilament number resulted in a predicted areal constriction ratio (the ratio of the internodal axonal cross-sectional area to the nodal axonal cross-sectional area) consistent with experimental observations on mouse sciatic nerve axons in vivo (Reles and Friede, 1991).

To test for consistency in our model, we compared the predicted outcomes of simulated pulse-escape experiments in the nodal and internodal domains with those observed experimentally in Walker et al. (2019). While not identical, we found that the predicted fluorescence decays in the node and internode were qualitatively similar, and within the error range of the experimental data. This is notable because the neurofilament kinetics in the node and internode in our model differ only by a single rate, the on-rate, which describes differences in neurofilament access to microtubule tracks. Thus, we conclude that the proximity of neurofilaments to microtubules is a potential regulator of neurofilament transport in axons and is sufficient to explain the local acceleration of neurofilaments through these axonal constrictions, ensuring a stable morphology across these physiologically important structures.

A notable feature of our model is that it predicts an interesting interdependency between the neurofilament flux j, microtubule number N_MT, neurofilament velocity , and axonal cross-sectional area A.The caliber is dependent on the neurofilament content, which is determined by the neurofilament influx from the cell body and the average neurofilament velocity. The average velocity is dependent on the ratio of the numbers of neurofilaments and microtubules, which is influenced by the neurofilament influx. Doubling the neurofilament influx and doubling the number of microtubules results in a doubling of the cross-sectional area without any change in the ratio of neurofilaments to microtubules and thus without any change in the average neurofilament velocity. However, if the number of microtubules is held constant, even a small increase in the neurofilament flux can have a comparably large effect on the axonal cross-sectional area. For example, for the internodal parameters and microtubule number used in this study, increasing the influx of neurofilaments by only about 12% from 7.8 to 8.7/s results in a decrease of the neurofilament transport velocity by 55% from 0.42 mm/d to 0.19 mm/d and a 2.5× increase in the number of neurofilaments from 213,500 to 533,750, leading to a doubling in the cross-sectional area from 10 to 20 μm². This underlines the critical and in general nonlinear influence of neurofilament flux on the axon cross-sectional area: a small change in flux can result in a large change in axon caliber.

We note that our model is based on the fundamental assumption that neurofilaments perform a diffusive search to bind to available microtubule tracks. This assumption that neurofilaments can diffuse within the radial dimension of axons is consistent with prior reports that neurofilaments behave as weakly interacting elements that distribute randomly in axonal cross-sections across a range of densities (Price et al., 1988) and diffuse apart from each other freely when separated from their plasma membrane (Brown and Lasek, 1993). However, this diffusion coefficient has not been measured experimentally. Therefore, the value of this parameter in our simulations was selected by matching the resulting neurofilament on-rate to within an order of magnitude of that predicted in previous work (Trivedi et al., 2007; Jung and Brown, 2009; Walker et al., 2019). Simulations with lower diffusion coefficients and all other rate constants unchanged led to smaller average numbers of on-track neurofilaments, as well as smaller on-rates and velocities, and larger nodal constrictions (unpublished data). Given the importance of the radial mobility of neurofilaments in our model, the development of methods that can measure the radial diffusion coefficient of neurofilaments in axons experimentally must be a priority for future experimentation.

Our model was designed to test the specific hypothesis at hand in a computationally efficient manner and therefore includes a number of simplifying assumptions that should be validated in future experimentation. One assumption in our model is that the neurofilament length distribution measured in cultured neurons applies to myelinated axons in vivo. We investigated the influence of this assumption in our simulations and found that the average neurofilament length influences the sharpness of the decline in neurofilament content. The longer the neurofilaments, the more gradual the decline. To test this prediction experimentally, it will be necessary to measure the neurofilament length distribution in myelinated axons in vivo and also perform fine-scale mapping of neurofilament number across nodes in these axons. A second assumption in our model is that the axon can be described as a one-dimensional domain along which transport processes occur, incorporating the neurofilament search for microtubule tracks in the radial dimension of the axon through the proposed diffusion model. This hybrid modeling strategy has the advantage of reduced computational cost and fewer unknown parameters. However, a more accurate approach would include the direct simulation of neurofilament motion within the crowded three-dimensional environment of the axon, incorporating the mechanical effects of the nodal morphology on neurofilaments navigating the constricted node. We plan to pursue this direction in future work, which will require experimental measurements of neurofilament polymer mechanics, organization, and interactions that are currently unavailable.

Another assumption in our model is that the axonal plasma membrane deforms freely to accommodate changes in neurofilament content without resistance. A more realistic model might include a viscoelastic boundary that exhibits elastic resistance on short time scales and viscous deformation on longer time scales. However, such a model would require measurements of the deformability of the axonal plasma membrane, which will be influenced by the mechanical properties of the membrane cytoskeleton, myelin sheath, and extracellular matrix. In practice, the assumption of a freely deformable boundary may be a reasonable approximation on the slow time course of axonal expansion and contraction caused by changes in neurofilament transport (i.e., hours or days). Moreover, this may not be of great consequence for our steady state simulations of nodes of Ranvier, where the fluctuations in neurofilament content are small. In future work, where we plan to simulate internodal axon expansion and nodal formation during development, a more detailed model of the membrane mechanics will be required.

It is important to reiterate that nodes of Ranvier and nodal constrictions are observed in mutant mice that lack axonal neurofilaments, though both the nodes and internodes in these mice fail to attain normal caliber (Perrot et al., 2007). This is consistent with the known role of neurofilaments as space-filling structures that expand axonal caliber, but it also indicates that neurofilaments are not required for nodal constrictions to form. Thus, we do not propose that nodal constrictions arise as a consequence of the local modulation of neurofilament transport, but rather that the local modulation of neurofilament transport is essential to allow neurofilaments to navigate these constrictions without piling up on either side. In fact, nodes are complex and highly structured domains with a distinct membrane architecture that are assembled independent of neurofilaments, triggered by interactions with the myelinating cells, and likely templated by a highly organized and periodic submembrane actin cytoskeleton (Susuki et al., 2016; Ghosh et al., 2018). How the expansive forces generated by neurofilaments interact mechanically with this nodal cytoarchitecture during axonal development is an intriguing question.

Neurofilaments are of clinical interest because they can accumulate excessively in a variety of toxic neuropathies and neurodegenerative diseases, leading to swelling and distortion of axons and consequent disruption of nerve conduction (De Vos et al., 2008; Millecamps and Julien, 2013). In fact, as we have noted previously (Walker et al., 2019), studies on animal models of neurodegenerative and neurotoxic disease have often reported that these accumulations appear proximal to nodes of Ranvier (e.g., Griffin et al., 1982; Jones and Cavanagh, 1983; Jacobs, 1984; Gold et al., 1986; Hirai et al., 1999; Lancaster et al., 2018). Our model suggests that this could arise due to a partial local destabilization of axonal microtubules resulting in a reduction in microtubule number, such as has been implicated in a number of neurodegenerative diseases. The requirement that neurofilaments accelerate locally in nodes to maintain a steady state morphology across these sites may make nodes particularly vulnerable to such perturbations.

FOOTNOTES

Abbreviations used:

inter	internode
MT	microtubule
NF	neurofilament
paGFP	photoactivatable green fluorescent protein
SD	standard deviation

ACKNOWLEDGMENTS

REFERENCES