A possible mechanism for neurofilament slowing down in myelinated axon: Phosphorylation-induced variation of NF kinetics

Neurofilaments(NFs) are the most abundant intermediate filaments that make up the inner volume of axon, with possible phosphorylation on their side arms, and their slow axonal transport by molecular motors along microtubule tracks in a “stop-and-go” manner with rapid, intermittent and bidirectional motion. The kinetics of NFs and morphology of axon are dramatically different between myelinate internode and unmyelinated node of Ranvier. The NFs in the node transport as 7.6 times faster as in the internode, and the distribution of NFs population in the internode is 7.6 folds as much as in the node of Ranvier. We hypothesize that the phosphorylation of NFs could reduce the on-track rate and slow down their transport velocity in the internode. By modifying the ‘6-state’ model with (a) an extra phosphorylation kinetics to each six state and (b) construction a new ‘8-state’ model in which NFs at off-track can be phosphorylated and have smaller on-track rate, our model and simulation demonstrate that the phosphorylation-induced decrease of on-track rate could slow down the NFs average velocity and increase the axonal caliber. The degree of phosphorylation may indicate the extent of velocity reduction. The Continuity equation used in our paper predicts that the ratio of NFs population is inverse proportional to the ratios of average velocity of NFs between node of Ranvier and internode. We speculate that the myelination of axon could increase the level of phosphorylation of NF side arms, and decrease the possibility of NFs to get on-track of microtubules, therefore slow down their transport velocity. In summary, our work provides a potential mechanism for understanding the phosphorylation kinetics of NFs in regulating their transport and morphology of axon in myelinated axons, and the different kinetics of NFs between node and internode.


The analytical and numerical result for '6-state' model with phosphorylation kinetics
2.1.1 The analytical calculation for the NF velocity with phosphorylation kinetics in each six states. We use mathematical analysis to predict the velocity of NFs at node and internode. When NFs transport reach dynamically equilibrium state, and it leads to: @P ph @t ¼ À g de P ph þ g ph P de ¼ 0 @P de @t ¼ þg de P ph À g ph P de ¼ 0 Therefore, the ratio of the probability of NFs at phosphorylated and dephosphorylated states is P ph : P de = γ de : γ ph .
The final average velocity of NFs is contributed by two parts: one is the velocity of NFs at the phosphorylated state, and the other is the velocity of NFs at the dephosphorylated state. Therefore, an equation can be written as the summation of the following two terms as in Eq (1): the first term is the product of the fraction of NFs at dephosphorylated state f de ¼ P de P ph þP de with the average velocity at the dephosphorylated state � v fast ; the second term is the product of the fraction of NFs at phosphorylated state f ph ¼ P ph P ph þP de , therefore the average velocity � v can be written as: In the following, we can calculate the average velocity in both node and internode sections along the axon. The probability ratio between dephosphorylation and phosphorylation is set as: P de : P ph ¼ g de : g ph ¼ 1: 8; x � 20mm and x � 30 mm; internode 8: 1; x 2 ½20; 30�mm; node of Ranvier ( ð1aÞ Then the average velocity at both node and internode can be calculated as: Where, � v fast ¼ 7:6� v; v node v internode ¼ 61:8 15:6 ¼ 3:96. This is the theoretic prediction. Our theoretical calculation successfully predicts the simulation result as shown in Fig 1. This result is based on the experimental results of � v fast ¼ 7:6� v as shown in [19], it can be generalized to other type of nerves or species according to our model in Eq (1). In addition, the rate of γ de : γ ph could be different for different nerves or species, therefore, the result can also be different according to Eq (1). This model is a general working frame that could be applied to different situations.

Monte
Carlo simulation with phosphorylation kinetics in each state of '6-state' model. To simulate the effects of phosphorylation on the transport dynamics of NFs and axonal morphology, we construct a single axon with 50 μm long and the section of node of Ranvier is from 20 μm and 30 μm. The length of node of Ranvier ranges around 1 or 2 μm, here we assume that the node length is 10μm in order to illustrate the simulation results clearly; and our model and simulation are independent of the length of node, but depend on the kinetics of NFs phosphorylation and other parameters related with transport, and therefore can be applied to any type of axon or species.
We add the kinetics of phosphorylation and dephosphorylation into each state of '6-states'. Based on Eq (6), we set the corresponding value of γ on in node and internode. At the initial section of the axon, we input NFs continuous into the axon (input one NF every 10s). We set the NFs at the initial of axon are at the state of running anterograde, if the displacement of the NFs is beyond 0-50 μm, we will remove it from our system. After reaching the steady state, the distribution of NFs population along the axon and their kinetic states are observed and recorded.
To improve the accuracy of simulation, the experiment was repeated 300 times. And we count the number of NFs in each position of the 300 axons, and calculate their mean value and the standard deviation (Fig 1).
In our simulation, we recorded each NF location and state information. Fig 1 shows the number of NFs at each position of the axon. We divided the 50 μm axon into 25 bins, where each bin is 2 μm long. We collected the number of NFs in each bin, and obtain the distribution of NFs in Fig 1. Besides, we can get the probability of NFs in each state along the axon.
According to continuity equation, flux ¼ N � � V, the flux of the continuously input NFs is calculated by the product of the number of NFs (N) with the average motion velocity ( � V). Therefore, when the NF transport reaches steady state. The equation of N node � � V node ¼ N internode � � V internode will be met. Based on this, we can test our simulation results. In order to test our simulation results, we can calculate the average velocity of NFs in the Monte Carlo simulation by the following formula: According to the results of Fig 1 and Eq (1b), we calculate the average velocity of NFs in the internode � V 1 ¼ 2:2 mm=day and in the Node of Ranvier � V 2 ¼ 8:7 mm=day and . This is simulation result, which is close to theoretical result of 3.96.
We can obtain the average number of NFs in internode section N 1 and in the node section N 2 by counting the number of NFs in each position along the axon in Fig 1, and the result is Remarkably, it can be seen that Therefore, we conclude that the continuity equation still holds true in this transport system, and we have the flux of NFs in the internode and nodes equals to each other: The velocity modulation depends on the rate ratio between phosphorylation and dephosphorylation r. The consistence of theoretical analysis and computational simulation has several indications.
The first indication is that, if the velocity ratio between node and internode has to be 7.6, then how much modulation of γ on has to be made? Given the relation and distribution of phosphorylation and dephosphorylation unchanged. Based on Eq (1), the reconstruction calculation can be done in the following.
We assume � v fast ¼ n� v, where n is the number unknown and needs to find out, according to Eq (1): If Eq (1a) stands, then it will lead to � V node : � V internode ¼ 8nþ1 nþ8 ¼ 7:6, subsequently this results in n = 149.5. This calculation indicates that if a velocity ratio between node and internode needs to be 7.6, the velocity of fast transport has to be boosted by 149.5 folds, i.e. � v fast ¼ 149:5 � v. According to construction of the on-track rates based on Eq (6) in the model and method section, the on-track rate can be calculated out.
In addition, the ratios of average velocity between internode and node can be written as functions of n, q 4 , where q 4 ¼ g de g ph , If it is set r ¼ g ph g de internode ð Þ, and assume in the node: g ph g de ¼ 1 r , then we have: According to Eq (1c), the ratio of velocity depends on n and r values in the node and internode. Therefore, if the two factors can be measured in experiments, then the average velocity between internode and node can also be predicted. We plot the velocity ratio between node and internode as function of n in Fig 2, where x axis is the velocity ratio n without the influence of phosphorylation kinetics, r is the ratio between phosphorylation rate and the dephosphorylation in the internode section. We assume in the node the ratio is inverse of it in the internode section, which is 1/r. For different r values, the colored curve show different dependence of y on the n and r, the velocity ratio between node and internode is increasing with r = γ ph : γ de for a fixed value of x. As is shown in the red curve (r = 8.0) in Fig 2, the x axis is only for 1 to100 but not reaches 149.5, but trend can be predicted there that for x = 149.5, the y axis can only reaches 7.6. In another word, if r = 9.0 in the dashed blue curve, the x axis would be x = 48.2 for the value of y = 7.6; if r = 10.0 in the dashed green curve, then x = 31.3 for the y = 7.6. Therefore, the ratio r needs to increase to reduce the effort of boosting velocity without phosphorylation kinetics.
The second indication is, we chose to increase γ on in the node to make the velocity faster than the internode, and the supposed boost of velocity is 7.6 fold; while the phosphorylation kinetics interfere the '6-state' model of NFs, by dividing each state into phosphorylated and dephosphorylated states, with different transition rates in the node and internode, the originally constructed velocity ratio between node and internode of 7.6 drops to be only 3.96. In another word, it demonstrates that the phosphorylation kinetics has decreased the influence of on-track rate γ on in determining the average velocity ratio of NFs between node and internode. The dependence of velocity boost on the rate ratio between phosphorylation and dephosphorylation r = γ ph : γ de . The x-axis is the velocity ratio n between the node V(fast) and internode V(slow) without the phosphorylation kinetics. The y-axis is the velocity ratio between node V(node) and internode V(internode) with the phosphorylation kinetics. https://doi.org/10.1371/journal.pone.0247656.g002

PLOS ONE
As can be seen in Fig 2, for all values of r, the curves are all going to be saturated at different values for extremely large x. It indicates that no matter how fast the velocity in the node is boosted compared to the internode, with a certain phosphorylation kinetics added into the system, the velocity ratios between node and internode will eventually approaches a saturated value, which is much smaller than the value n without phosphorylation kinetics. The higher values of r, the higher the saturated ratio is going to be. In general, for all values of r, the y axis is still below the y = x curve;, which indicates that the phosphorylation kinetics has dramatically decreased the velocity ratios between node and internode with the phosphorylation kinetics added into the system.
The last point is that, as the ratio r changes, the velocity ratio between the node and internode will be regulated. We take an extreme example to illustrate the effect. If r = γ ph : γ de = 1 in the internode, then the average velocity will become � On the other hand, if γ ph � γ de in the node, then the average velocity will become In summary, we can see that the at the limit of extremely phosphorylation in the internode and no phosphorylation in the node, the velocity ratio will become to be � In the other extreme, if the phosphorylation ratios r = γ ph : γ de = 1 in the internode, then the average velocity will also be the same between node and internode � V internode : � V node ¼ 1: 1.

Simulation result of '8-state' model 2.2.1 Simulation result of reconstruction of velocity in '8-state' model by numerical simulation of PDE.
The partial differential equations in Eq (8) are discretized into ordinary differential equation by First Order Upwind to obtain numerical solution. The on-track rates reconstructed in Fig 8 are put into the simulation. In the simulation, NFs are injected into the proximal section of axon with total 50 μm long. The boundary conditions are set as: P� (x 0 , t) = 0; P� (x l , t) = P� (x l − 1, t); x 0 and x l denotes the proximal boundary axon and the distal boundary of axon, P� denotes each probability of eight states. Initial condition of the PDE of NFs distribution is set by the equilibrium distribution according to Eq (9).
In order to observe the distribution of NFs along the axon in the equilibrium state, NFs are continuously input into the axons from the proximal section (input one unit of NF population every 1s). The results of the numerical simulation are shown in the Fig 3. As can be seen, the NFs distribution along axon does not change after 9 hours, reaching equilibrium state. Note that it is a dynamical equilibrium, with constant flux of NFs entering into the proximal axon and the same constant flux exiting axon at the distal section.
As can be seen in Fig 3, the ratio of population of NFs between internode and node of Ranvier is 2.415:0.3175 = 7.6. Note that the section form 20 μm to 30 μm is the node of Ranvier. It not only reproduces experimental observation results [6], but also consists with the theoretical prediction in Sec.4.2.1.
As can be seen from Fig 4, the probability of NFs at each eight state is differentially distributed along the axon. When the NFs enter the node of Ranvier, due to the higher dephosphorylation rate and thus higher on-track rate, probabilities of NFs at on-track states P a , P a0 , P r , P r0 increase, while the probability of P ap , P rp decrease. At the node of Ranvier, the NFs at off-track states are dephosphorylated, thus NFs at P ap1 , P rp1 are much lower; while NFs in the internode are highly phosphorylated in the axon, thus the NFs at states of P ap1 , P rp1 are much higher in the internode. Therefore, the introduction of transition kinetics with different rates of phosphorylation and dephosphorylation in node and internode, and its modulation on the "ontrack" rate g 2 on can eventually demonstrate the influence of phosphorylation on NF transport.  According to the results of Fig 5 and Eq (1b), we calculate the average velocity of NFs in the internode, From the simulation result of Fig 5, we can count the average number of NFs in the internode N 1 = 10.82 and node of Ranvier N 2 = 1.44 and we have N 1 N 2 ¼ 10:48 1:37 ¼ 7:51, which is close to the inverse ratio of average velocity: Therefore, the reconstruction is valid by holding the continuity equation true along the same axon of node and internode sections.
In addition, this reconstruction fully considered the phosphorylation induced on-track rate change, and the simulation results from PDE and Monte Carlo method are the same. Our method successfully reproduces the experimentally observed NFs population distribution along axon. Our model provides a possible mechanism for understanding the NFs slowing down in the myelinated internode.

Analysis of the on-track rate g 2
on and phosphorylation rate ratio r. If the rate ratio between phosphorylation and dephosphorylation in the internode r = γ ph : γ de (and the ration will be 1/r at node of Ranvier) varies, and the on-track rate g 2 on changes simultaneously for both internode and node, while the rate g 1 on ¼ 0:000275s À 1 is unchanged, then it can be seen how the velocity ratio between node and internode changes with the parameter r in Fig 6. The larger r value is, the higher velocity boost in the node will be for a certain on-track rate g 2 on . For one certain value of r, there is always an optimized value of g 2 on to maximize the velocity ratio in the y axis. In order to reach a desired value in y axis (boost the velocity in the node), it may not work by simply increasing g 2 on , the value of g 2 on has to cooperate together with parameter r to reach the desired velocity acceleration in the node. Another point from Fig 6 is that, if g 2 on is the same for node and internode, then it is not possible to reach the desired velocity ratio in the y axis for a chosen r. For example, if r = 8 in the red solid curve, and no matter how much variation of g 2 on has, it is not possible to reach a velocity ratio 7.6. In Fig 7, for the case of same on-track values g 2 on for both node and internode, as is shown in the blue curve (gon2_node = b � 1), the value of r has to be more than 30.0. to reach the velocity ratio of 7.6 in y axis.
On the other hand, as the ratio of on-track rate g 2 on between node and internode increases, the velocity ratio between node and internode can be boosted for different values of r, as can be seen in Fig 7. Another point can be seen in Fig 7 is that, when the g 2 on in the node is smaller than the internode, the velocity ratio can also be larger than one if the r value is selected poroperly. For example, for the dashed red curve, the g 2 on value in the node is half of it in the internode, the y axis can still be larger than 1 and can reach 5 for larger r values. This indicates that the phosphorylation kinetics itself can balance the lower on-track rate in the node and still result in a faster velocity in the node of Ranvier.

Conclusion
Based on the "stop-and-go" hypothesis, firstly we tried to add phosphorylation kinetics in each state, and analyzed the effect of phosphorylation on the average velocity of NFs transport in the node and internode. The basic assumption is the on-track rate is modulated by the phosphorylation kinetics, and NF transport velocity and the distribution of NFs population are thus changed by the on-track rate, according to the continuity equation. Secondly, by assuming that only at off-track states, NFs have phosphorylation kinetics and the phosphorylation status will eventually regulate the on-track rate, in which we introduce the kinetic transition between phosphorylation and dephosphorylation of NFs and build the '8-state' model. Our model demonstrated that the regulation of phosphorylation on NF transport and axon morphology, and provided a potential mechanism for NFs slowing down during myelination.
The modulation of the on-track rates of NFs is based on the assumption that NFs in node of Ranvier are less phosphorylated thus the on-track rate is higher, while in the internode, due to highly phosphorylation, NFs on-track rate is lower. Another basis of our modulation is that the NFs flux is a constant at equilibrium state, which results in the inverse proportional relationship between population of NFs and the average velocity distribution along node and internode. Through the Monte Carlo simulation and numerical solution of PDE of '8-state' model, our results demonstrate the dephosphorylation of NFs can accelerate NFs kinetics in the node of Ranvier and phosphorylation of NFs can decrease NFs kinetics in the internode. In the end, our simulation results show that ratio of the number of NFs is inverse proportional to the average of velocity of NFs at both node and internode. Our model provides a potential mechanism to understand the regulation of NFs kinetics at myelinated axon by NFs phosphorylation kinetics; besides, the dynamic equilibrium for the kinetics between phosphorylation and dephosphorylation determines the fraction of NFs at each state, which regulates the ontrack rate of NF transport, and this mechanism can be applied to any axon for any species.

Discussion
However, our model still needs to be improved in the following aspects. First, we assume that the NFs are independent of each other, and do not consider the interaction between them, which is certainly unrealistic in the real axonal transport. Besides, NF-NF association results in the formation of NF "bundles", the "bundles" may have important effect on the axonal transport. Second, the transition rate constants γ 10 , γ 01 , γ ra , γ ar are put into our model according to the "stop-and-go" model [18] and by the observation data of rat superior cervical ganglion neurons published by the research [19]. Noting that the node of Ranvier is relatively short compared with the internode section of axon, it must be difficult to measure the average velocity of NFs passing through the node in the experiments. When the experiment has new progress, our model can modify these parameters in according to the new experimental data.
Our model assumed different rates of phosphorylation and dephosphorylation between nodes and internodes, which have not been experimentally observed, and may have different values in different species or types of nerves. However, Data presented in [38] clearly distinguish between the densities of NF phosphorylated epitopes in internodes versus nodes of Ranvier. NF-H and NF-M phosphorylated epitopes are reduced by 60 and 40% respectively in the nodes relative to internodes. The basic principle of our model can predict transport velocities of NFs based on different values of phosphorylation rate constants. In addition, according to the literature, the phosphorylation rates of NFs may change under differential myelination level [24,27,39]. In this sense our model provides a possible mechanism for understanding myelination regulation of NFs transport inside the axon.
There are different results, regarding on the experimental observation of NFs transport due to the NFs phosphorylation. Some experiments show that the phosphorylation slow down the NFs transport [1,40,41], while others show that phosphorylation did not slow down NFs transport [34][35][36][37]. We can interpret that the level of phosphorylation and dephosphorylation in regulating the NFs on-track rate might be the key. For example, if the NFs in the internode are heavily phosphorylated, it might decrease the on-track rate severely, then the average velocity will be slowed down. On the other hand, if NFs are less phosphorylated in the internode, the on-track rate will be less affected, and the average velocity will not be slowed down. During the myelination process occurring in the developmental period, the caliber of axon gradually increases and it might be interpreted as due to the gradual phosphorylation of the NFs which decreases the on-track rates gradually, and results in the slowing down of NF transport. Therefore the phosphorylation effect can be both spatial and temporal, and the effect on the rate of on-track and NFs average transport velocity will follow, and that needs to be investigated in the future work.

Phosphorylation kinetics added to each state of '6-state' model
In this section, we add NF transition kinetics of phosphorylation and dephosphorylation into each state in the '6-state' model. By assuming that phosphorylated NFs have lower on-track rate and dephosphorylated NFs have higher on-track rate, we show that NFs phosphorylation regulates the average velocity of the NFs and their final distribution along axon by both analytical solution and simulation results.

Model description of the phosphorylation kinetics added into the '6-state' model.
We model NFs transport based on a previous model: '6-state' model [18], where each NF moves bi-directionally along the axon switching between six kinetic states (P a , P r , P a0 P r0 , P ap , P rp ). The model gave an analytical expression for average velocity expression; and the model can produce simulation results in a good approximation of experimental data for the group transport wave of Gaussian function. The model mathematical expression can be written by coupled Partial Differential Equations (PDE) in Eq (2). @P a @t ¼ À v a @P a @x À g 10 P a þ g 01 P a0 @P r @t ¼ À v r @P r @x À g 10 P r þ g 01 P r0 @P a0 @t ¼ À ðg 01 þ g ar ÞP a0 þ g 10 P a þ g ra P r0 þ g on P ap À g off P a0 @P r0 @t ¼ À ðg 01 þ g ra ÞP r0 þ g 10 P r þ g ar P a0 þ g on P rp À g off P r0 @P ap @t ¼ g off P a0 À g on P ap À g ar P ap þ g ra P rp @P rp @t ¼ g off P r0 À g on P rp À g ra P rp þ g ar P ap Where P a (x, t), P r (x, t), P a0 (x, t), P r0 (x, t), P ap (x, t), P rp (x, t) represent the distribution of running-anterograde, running-retrograde, pausing-anterograde, pausing-retrograde, off-trackanterograde, and off-track-retrograde, respectively. In our simulation, the anterograde and retrograde velocities v a = 0.52 μm/s, v r = −0.36μm/s are adopted from research [2,7,8,18,40,41]. The spatial distribution of NFs population is: If NFs enter the axon at the proximal section and leave the axon at distal section at a constant rate, after long time, the steady distribution of NFs in the kinetic states can be obtained readily by solving the above equations. The theoretical solution of probability at each of the six states are [18]: For the paper, the rate parameters are set as [18]: γ 10 = 0.14 s −1 , γ 01 = 0.064 s −1 , γ off = 0.00445 s −1 , γ ra = 1.4 × 10 −5 s −1 , γ ar = 4.2 × 10 −6 s −1 . The rate of on-track will be modified in the latter section.
According to the experimental result, the average movement velocity of NFs in the node of Ranvier is 7.6 times as that in the internode sections [19]. Why the NFs can accelerate so rapidly at node of Ranvier? Here, we introduce the transition kinetics between phosphorylation and dephosphorylation into the '6-state' model, and investigate how phosphorylation of NFs could regulate the kinetics of NFs. At first, we assume for each six state, NFs can have the kinetics of phosphorylation and dephosphorylation as seen in Fig 8. As NFs can have different episodes for phosphorylation either in NFH or in NFM or both [42,43], and the NFs phosphorylation at different episodes and how many episodes are phosphorylated involves complicated chemical reactions of kinase or other proteins [44], in order to make our model concise and can be generalized to various conditions, we from biophysical perspective assume that the NFs as a one unit that can have two state of phosphorylation, and the percentage of phosphorylated or level of phosphorylation can be recognized as the fraction of NFs at phosphorylated state compared to the dephosphorylated state.
The kinetics of phosphorylation and dephosphorylation can be written in the following equation: @P ph @t ¼ À g de P ph þ g ph P de @P de @t ¼ À g ph P de þ g de P ph Where P ph (x, t) and P de (x, t) represent the distribution of phosphorylated and dephosphorylated NFs respectively; the rate of γ de and γ ph are assumed to be dependent on the location where NFs transport through, which is related to myelination function along the axon. According to the experimental observations [18,25,30,32], we assume that the NFs has a higher rate of phosphorylation at internode due to the myelination process around axon, while NFs has a higher rate of dephosphorylation at node of Ranvier, where there is no myelination wrapped, as can be seen in Fig 9. The rate of phosphorylation and dephosphorylation in the internode is γ ph = 0.8, γ de = 0.1. while the node, the rates are set opposite: γ ph = 0.1, γ de = 0.8.

Reconstruction of the on-track rate.
We assume that phosphorylation of NFs makes NFs difficult to get on-track of microtubule by either tangling with neighboring NFs or affecting the interaction between the NFs and molecular motor kinesin or dynein, therefore, therefore, it decrease the rate of on-track γ on . On the contrary, NFs dephosphorylation can reduce side arm length, and facilitate NFs to jump on-track for running, which results in higher on-track rate. In the following part, the on-track rate reconstruction is illustrated based on the continuity equation [14].

PLOS ONE
From paper [18], the average movement velocity � v is: The magnitude of � v depends on the coefficients q 1 , q 2 and q 3 . The on-track rate of γ on can be re-written as a function of average velocity as follows: Based on Eq (5), with other transition parameters unchanged, if the NF average velocity in the node of Ranvier is 7.6 times as that in the internode section, then the ontrack rate will become g 0 on ¼ 5:16 � 10 À 3 s À 1 (corresponding to average velocity of � v fast ¼ 7:6� v ¼ 7:6 � 0:71mm=day ¼ 5:31mm=day), while average velocity in the internode, the on-track rate is only g on ¼ 2:75 � 10 À 4 s À 1 ð� v ¼ 0:71mm=dayÞ.

PLOS ONE
Thus, we set the value of γ on to be location dependent as follows: The rate parameter about phosphorylation of NFs of γ de , γ ph along the axon is shown in Fig  9. At the node of Ranvier (20μm < x < 30μm), the rate of phosphorylation (γ ph = 0.1 s −1 ) is much lower than the rate of dephosphorylation (γ de = 0.8 s −1 ), indicating that NFs are less phosphorylated; while in the internode section (x � 20 μm and x > = 30 μm), the phosphorylation rate (γ ph = 0.8 s −1 ) is much higher than the rate of dephosphorylation (γ de = 0.1 s −1 ), indicating that NFs are more phosphorylated. In the constructed configuration, the model can demonstrate the effect of phosphorylation on the transport of NFs.

The newly developed '8-state' model
In this section, we build a new model-"8-state" model by considering two extra states of phosphorylated and dephosphorylated in the "6-state" model, and derive a mathematical expression for the average velocity of NFs.

Model description of '8-state' model.
Based on the experimental observations of axonal morphology and the NFs population difference between the internode and node of Ranvier [19], we assume that NFs transport decrease in the internode because of myelination and phosphorylation, while NFs in the node of Ranvier accelerate their transport, and our

PLOS ONE
model increases the on-track rate from the less phosphorylated state to on-track pausing state (g 2 on ) based on the continuity equation. Our model can successfully reproduce higher NFs content in the internode and lower distribution in the node of Ranvier by 7.6 folds. In addition, by the Monte Carlo simulation and PDE solution it is true that the continuity equation holds true for the '8-state' model.
Considering that NFs have transition kinetics of phosphorylation and dephosphorylation, we modify the '6-state' model by subdividing two off-track pausing states at both anterograde and retrograde directions, resulting in an ' 8-state' model.
The basic assumptions of our model is that we assume only NFs at off-track states have to be divided into two sub-states, (which is different from above, where the NFs at each of the six states have both phosphorylated and dephosphorylated states): one is NFs of phosphorylation with P 1 ap and P 1 rp ; the other is NFs of de phosphorylated with P 2 ap and P 2 rp ; and their on-track rate will be different due to the phosphorylated status, with lower value of on-track rate (g 1 on ) in phosphorylated sub-states and higher on-track rate in dephosphorylated state (g 2 on ). The schematic diagram of the transformation between the eight states is in Fig 10. Correspondingly, the kinetic equations for NFs transport can be written as follows: @P a @t ¼ À v a @P a @x À g 10 P a þ g 01 P a0 @P r @t ¼ À v r @P r @x À g 10 P r þ g 01 P r0 @P a0 @t ¼ À ðg 01 þ g ar þ 2g off ÞP a0 þ g 10 P a þ g ra P r0 þ g 1 on P 1 ap þ g 2 on P 2 ap @P r0 @t ¼ À ðg 01 þ g ra þ 2g off ÞP r0 þ g 10 P r þ g ar P a0 þ g 1 on P 1 rp þ g 2 on P 2 ap @P 1 ap @t ¼ g off P a0 À ðg 1 on þ g ar þ g de ÞP 1 ap þ g ra P 1 rp þ g ph P 2 ap @P 2 ap @t ¼ g off P a0 À ðg 2 on þ g ar þ g ph ÞP 2 ap þ g ra P 2 rp þ g de P 1 ap @P 1 rp @t ¼ g off P r0 À ðg 1 on þ g ra þ g de ÞP 1 rp þ g ar P 1 ap þ g ph P 2 rp @P 2 rp @t ¼ g off P r0 À ðg 2 on þ g ra þ g ph ÞP 2 rp þ g ar P 2 The spatial distribution of NF population is the summation of all the motion states of individual NF, which can be written as Pðx; tÞ ¼ P a ðx; tÞ þ P r ðx; tÞ þ P a0 ðx; tÞ þ P r0 ðx; tÞ þ P 1 ap ðx; tÞ þ P 2 ap ðx; tÞ þ P 1 rp ðx; tÞ þ P 2 rp ðx; tÞ If NFs enter the axon at the proximal end and exit the axon at a constant rate, the equilibrium distribution of NFs can be reached, and can be obtained readily by solving the above Eq (8). The Eq (8) approximates Gaussian distribution. After derivation, we can obtain the similar formula for the average velocity and steady state distributions as in Eq (2), but with a different q 2 : where a 1 ¼ g ph g ph þ g de ; a 2 ¼ g de g ph þ g de P ap ¼ P 1 ap þ P 2 ap ; where; P 1 ap ¼ a 1 P ap ; P 2 ap ¼ a 2 P ap P rp ¼ P 1 rp þ P 2 rp ; where; P 1 rp ¼ a 1 P rp ; P 2 rp ¼ a 2 P rp 4.2.2 Velocity reconstruction for the node and internode in '8-state' model. The reproduction of NFs distribution in node and internode can readily adopt the strategy of Continuity Equation [14], since we have already derived the analytical expression of average velocity � v in Eq (9). It can be exemplified by reconstruction of the on-track rate in node and internode respectively, and then numerically solving the partial differential Eq (8) to reproduce the distribution of NFs along axon for both node and internode.
The on-track rate g 2 on at node and internode are reconstructed by the average velocity Eq (9), by keeping the on track rate g 1 on unchanged g 1 on ¼ 2:75 � 10 À 4 s À 1 , then we have: For the internode, the average velocity is assumed to be � v internode ¼ 1:01 mm=day, therefore, in the internode the rate parameters are: g 2 on ¼ 5:16 � 10 À 3 s À 1 ; the phosphorylation rates are set as before γ ph = 0.8 s −1 , γ de = 0.1 s −1 .
Therefore, the average velocity of NFs at node and internode is the same as we constructed in the beginning: � v node � v internode ¼ 0:0890mm=s 0:0117mm=s ¼ 7:6. The PDE simulation result is the same as the analytical calculations, as can be seen in the Fig 3. Supporting information S1 Data. (RAR)