A nervous system is composed of interconnected neurons, each of which has multiple short dendrites that receive signals from upstream neurons and a single long axon that transmits signals to downstream neurons. Early in development, however, a neuron has multiple short neurites of similar lengths that extend and retract repeatedly and apparently randomly (coined by [1] as Stage 2). These neurites later differentiate into dendrites and a single axon (Fig. 1). The symmetry-breaking process by which a single axon emerges from the neurites (Stage 3) is referred to as neuronal polarization. Surprisingly, this process does not require release of growth factors from target cells, as neuronal polarization has been shown to occur within in vitro hippocampal neuron cell cultures in which there are no growth factor gradients guiding the selection process [1,2]. The symmetry breaking appears to be random, and in experiments where the initial axon was cleaved off, a new one emerged from one of the other neurites [3]. Subsequent experimental manipulations demonstrated that it was possible for any neurite to become the winning neurite in the winner-takes-all contest [4,5]. It has been shown that even in the in vivo setting where growth factor gradients are present, the developing neurons go through these stages, lasting about a day, characterized by growth and retraction of neurites with a subsequent symmetry breaking event [6,7]. In this case, growth factor gradients influence the selection process, and indeed the random growth and retraction of neurites is thought to be a way for the neurites to explore the environment to seek out the growth factors (called neurotrophins). Nonetheless, even in vivo there is a winner-takes-all process that takes place, with bias provided by neurotrophin gradients.
There have been many studies aimed at understanding the biophysical mechanism of the neuronal polarization process [4,[8][9][10][11][12][13]. Several potential mechanisms have been identified, and it is clear that the process involves positive feedback signals to promote the growth of the axon as well as negative feedback signals to prevent the emergence of a second axon [10,[13][14][15][16]. In this article, we focus on the polarization process that occurs in vitro, without external neurotrophin gradients, using a minimal model for the winner-take-all selection process that incorporates positive feedback and explore the efficacy of several different negative feedback mechanism in generating persistent neuronal polarization.
The model is constructed based on the hypothesis that the dynamics underlying neurite growth and retraction are such that the system is bistable. That is, each neurite has two stable equilibria, short and long. The bistability is a product of the positive feedback. A second hypothesis is that the selection process is truly random, so no neurite is biased towards winning the competition to become an axon. We then explore three mechanisms of negative feedback.
One of these involves the retraction rate that is common to all neurites. The others involve a stochastic term that reflects randomly-timed and uniformly distributed actin waves which are known to be key to neurite elongation; each actin wave provides growth spurts by locally increasing the neurite volume to allow for microtubule polymerization [12,17,18]. We consider the effects of making neurite retraction, actin wave magnitude, and actin wave frequency dependent on the combined length of the neurites such that increased length increases the retraction rate, or decreases the actin wave magnitude or frequency. In all cases, the negative feedback is unbiased.
The results demonstrate that targeting the negative feedback to the stochastic growth magnitude (i.e., the actin wave term) results in the most persistent polarized system. They also demonstrate that having more than 2 neurites, but less than some upper bound, is optimal for achieving and maintaining neuronal polarization. This is consistent with the finding that most nascent neurons have between 2 and 10 neurites [2]. One model for neuronal polarization is based on competition for limited resources, including growth factor [16]. In the last section of Results, we demonstrate that a simplified model based on this limited-resource model contains the two elements that we find to be most successful at achieving persistent polarization: bistability and length-dependent reduction in the amplitude of actin-wave-driven stochastic excitation .
We consider a small population of R neurites, each with length L i , i = 1, 2, . . . , R. The basic model contains a term for positive feedback, a retraction term, and a stochastic term reflecting randomly-timed actin waves. The negative feedback is included later. The basic model is:
The first term reflects positive feedback through intracellular signaling [8,10,11,13,15,19], length-dependent diffusion of polarity effectors [16,20,21], and stablization of microtubules [22,23]. Positive feedback is an essential ingredient of bistability [24]. The second term provides a constant rate of neurite retraction, reflecting the retraction that occurs in all neurites between the arrival of actin waves [25]. The last term includes a sum of delta functions that describes sudden neurite elongation due to actin waves [25][26][27]. Each wave induces a jump in length of size A. The term t
is the time when the tip of the ith neurite receives the nth wave, which follows a Poisson process of rate λ (this is also the average number of waves generated per unit time).
With our minimal model, a neurite can be thought of as a particle in a double-well potential, as shown in Fig. 2A.
The left potential well corresponds to the state of being a short neurite, and the right well corresponds to the state of being a long neurite that will become an axon. The actin waves then provide random excitations that can drive the particle across the potential barrier at location L b , marking the establishment of an axon. All nascent neurons have short neurites, so they begin in the left potential well.
Since only one neurite develops into the axon of a typical neuron, the growth of other neurites should be suppressed to avoid having multiple axons. One way to incorporate this negative feedback into the model is to have the retraction rate r increase as the neurites get longer. For simplicity, we use the following length-dependent retraction rate:
where r 0 is a basal retraction rate and α is a parameter that controls the degree of suppression. The retraction rate is the same for every neurite, so the suppression is unbiased. Reducing the growth rate g in an unbiased way will give qualitatively similar results. Targeting the negative feedback to the retraction rate has the effect of eliminating the upper equilibrium state for each of the short neurites (Fig. 2B). Biologically, the increased retraction rate reflects the collection of inhibitory signals in the cell body sent from the neurite tips.
Since actin waves drive neurite growth, negative feedback can also be implemented by suppressing the generation of the waves. In our model, this is done by reducing the excitation rate λ according to the following equation:
where λ 0 is a basal excitation rate and µ controls the degree of reduction. Again, the inhibition is unbiased because the actin waves are shared equally among neurites. Under the rate reduction, a neurite retracts significantly between two waves, so its net growth is small (Fig. 2C). A length-dependent decrease in actin wave frequency is consistent with the observation that actin waves are less frequent once an axon is formed [18]. This inhibitory mechanism was also implemented in a previous mathematical model [21,28]. Finally, we implement negative feedback by reducing the amplitude A as follows:
where A 0 is a basal excitation amplitude and ϕ controls the degree of reduction. In terms of the particle in a double well potential, a reduced amplitude means that more excitations will be required to cross the potential barrier (Fig. 2D).
Biologically, actin waves carry growth factors produced at the cell body, so amplitude reduction could reflect depletion of the growth factors. A similar amplitude reduction mechanism was adopted in a previous modeling study [16].
For mathematical simplicity, we consider a nascent neuron with two neurites (R = 2) in this section. In fact, neurons with two neurites were also observed in experiments [2]. For such a neuron, we study the effect of each of the three negative feedback mechanisms from three perspectives: (1) the joint probability density of the lengths, denoted by p(L 1 , L 2 ), (2) the underlying deterministic phase portraits, and (3) the stochastic dynamics. To analyze the probability density and stochastic dynamics, we employ two complementary methods: the generalized cell-mapping method (GCM) and Monte Carlo (MC) simulations, which are explained in detail in the Appendix. The GCM allows us to efficiently determine the probability density of the lengths and its long-term limit. However, when dealing with neurons possessing more than two neurites, the GCM becomes computationally expensive. In such cases, the MC method proves to be more efficient, particularly when the timescale is short. Additionally, the MC method unveils neurite dynamics that are not captured by the probability density obtained through the GCM. Nevertheless, the MC method is less effective than the GCM in analyzing the long-term behavior of the probability density.
We begin by considering negative feedback through a length-dependent increase in the neurite retraction rate (Eq. ( 2)). The long-term joint probability distribution of the lengths exhibits two peaks, which indicates that the system spends most of the time near these peaks (Fig. 3A). Each peak represents a state with a long neurite and a short neurite, which we refer to as a polarized state. The formation of the peaks can be inferred from the following deterministic system:
where the term 1 2 Aλ is the time average of the Poissonian term in Eq. ( 1). The factor 1/2 accounts for the fact that the total number of actin waves is divided between the two neurites. We continue to use L i for the neurite lengths in this system, though they are no longer random variables. Figure 3B shows nullclines and the vector field for Eq. ( 5). There are seven equilibria, three of which are stable. Two correspond to the polarized state and are near the peaks of the probability distribution, while one near the origin corresponds to a state in which both neurites are short. The basin of attraction of the latter is small, so in the stochastic system described next trajectories leave this region quickly. To study the effect of the increased retraction alone, we set λ = λ 0 and A = A 0 . Other parameter values are given in Table 1.
The bimodal probability density does not necessarily imply a firmly established axon. Monte Carlo simulations show that both neurite lengths alternate between two levels. That is, the system makes frequent transitions between the two polarized states (Fig. 3C). This is further indicated by the short mean transition time from one peak to the other (O(10 2 ) time units; see Fig. 3A). Thus, a single polarized state is not maintained when the negative feedback is implemented upon the retraction rate.
We next explore the dynamics in which negative feedback is implemented through length-dependent reduction in the excitation rate as prescribed by Eq. ( 3). The joint probability density initially has two peaks at the two polarized states. As time progresses, however, the polarized peaks fade and a third peak corresponding to two long neurites gains prominence (the top right peak in Fig. 4A). The mean transition time from a polarized state to the nonpolarized state is O(10 3 ) time units (Fig. 4A). All three states appear as stable equilibria in the deterministic system, in addition to the stable equilibrium with a small basin of attraction corresponding to two short neurites (Fig. 4B). A Monte Carlo simulation shows the early development of a polarized state, followed by a transition to a state with two long neurites at t ≈ 1000 (Fig. 4C). These results indicate that this form of negative feedback is not effective at maintaining a persistent polarized state. Incorporating both length-dependent increased retraction rate and reduced excitation rate eliminates the two long-neurite state (Fig. 5A andB), but does not prevent flipping between polarized states (Fig. 5C).
Finally, we consider negative feedback implemented through a length-dependent reduction in the excitation amplitude as prescribed by Eq. ( 4). Unlike the excitation rate reduction, the amplitude reduction yields a bimodal probability density that develops almost immediately (at t = O(10), see Fig. 6A) and persists even at t = 10 6 (to be explained in the next section). A peak in which both neurites are long does not appear until much later, and the mean transition time from a polarized state to this nonpolarized state is O(10 9 ) (Fig. 6B), which is much larger than the mean transition time when negative feedback is through rate reduction. Biologically, this means that the polarized state persists long enough that later stages of neuron development, including targeting of the nascent axon to appropirate targets via neurotrophins, can occur. Also, the system does not flip between the two polarized states, as shown with a Monte Carlo simulation (Fig. 6C). These results indicate that implementing negative feedback through a length-dependent reduction in the excitation amplitude results in persistent neuronal polarization. Other parameter values are given in Table 1.
Figure 5: The dynamics of the two-neurite system with negative feedback implemented through both a lengthdependent increase in the retraction rate and a decrease in the excitation rate. (A) The long-term probability density is bimodal with two polarized states. The mean transition time between these states is O(10 3 ). (B) The phase portrait shows two stable polarized equilibria and a stable equilibrium with small basin of attraction in which both neurites are short. (C) A Monte Carlo simulation shows the system flipping between the two polarized states. We set α = 0.02 and A = A 0 here. Other parameter values are given in Table 1.
We demonstrated above that regardless of the target of the negative feedback, the system enters a polarized state for some time before exiting to either (1) a different polarized state (i.e., flipping) or (2) a nonpolarized state (i.e., a state in which both neurites are long). However, the time that the system is in the polarized state varies greatly with the different forms of negative feedback. In this section, we examine why the persistence of the metastable polarized state is so different with the different negative feedback mechanisms. For this, we employ a tool called the ϵ-committor, developed by Lindner et al. [29]. It provides an estimate of the probability that a stochastic trajectory remains in a region of phase space for a duration of 1/ϵ, where ϵ is the rate at which the trajectory is moved into an absorbing state connected to the region. A definition and description of the calculation of the ϵ-committor is given in the Appendix.
We focus on a region that encloses the upper left peak in any of the bimodal probability densities in the previous Figure 7A shows that when the negative feedback is on the excitation amplitude, Other parameter values are given in Table 1.
C ϵ ≈ 1 over a timescale of 10 8 , meaning that with high probability a trajectory starting in R S remains in R S during this period of time. In contrast, C ϵ falls to zero much earlier when the negative feedback is on the excitation rate.
This indicates that the system leaves R S quickly and explains the rapid transition from a bimodal probability density to a unimodal probability density corresponding to two long neurites (Fig. 4). C ϵ starts to drop even earlier when the negative feedback is on the retraction rate, and it reaches ≈ 0.4. This means that the system spends about 40% of time in R S in the long run, consistent with its flipping behavior. These ϵ-committor results demonstrate again that applying the length-dependent negative feedback to the excitation amplitude works best in maintaining a unique polarized state. The difference in the variations of C ϵ when negative feedback is on excitation amplitude versus rate can also be quantified via the mean escape time. That is, the time at which a trajectory in a polarized state escapes to the other polarized state or to the nonpolarized state. Consider a single neurite of length L that follows Eq. ( 1). We use the general cell-mapping method (see Appendix) to calculate the mean time that L, starting at L = 0, exceeds L b for different values of A and λ (mimicking the effects of negative feedback on either of these two targets).
We find that the mean escape time (denoted by T c ) increases faster as the excitation amplitude is reduced than when rate is reduced (Fig. 7B). Therefore, reducing the excitation amplitude is more effective than reducing excitation rate on keeping a trajectory within an attracting basin. This explains the long persistence of the polarized state with negative feedback upon excitation amplitude.
Finally, we estimate the probability of crossing the threshold L b starting from L = 0 for a single neurite. To overcome retraction, the neurite must receive at least L b /A excitations during a short period (for simplicity, we assume that L b /A is an integer here, which is true for the parameter values listed in Table 1. If L b /A is not an integer, we need to round it up to the nearest integer. But this won't affect our result qualitatively). Consider τ = 1/r, the timescale of retraction. Let P c be the probability of having L b /A excitations during τ , which follows a Poisson distribution:
With the Stirling's Approximation for factorial
we get
If the amplitude A is reduced to A/m (m > 1) and λ remains unchanged, then
Thus, ln(P c ) decreases faster than linear reduction. To see this more clearly, we plug in the parameter values in Table 1, namely
On the other hand, if A is unchanged and λ is reduced to λ/m (m > 1), then
With the parameter values mentioned above, we get
where we can see that ln(P c ) decreases with m logarithmically. Thus, the probability of crossing the threshold is much larger at small excitation rate compared with the probability at a small excitation amplitude. One can conclude from this that the persistence in a single polarized state is greater when negative feedback is applied to excitation amplitude. 1.
Most developing neurons have between 2 and 10 neurites [2]. We next focus on the most effective form of negative feedback, targeted to the excitation amplitude, in model systems with more than 2 neurites. Results are obtained through Monte Carlo simulations, since the GCM approach to obtaining probability distributions is computationally expensive at higher dimensions. We seek to determine how the number of neurites R, as well as excitation amplitude and frequency, impact the probability of obtaining a single persistent polarized state.
The first set of results shows that the probability of obtaining a persistent polarized state first rises and then falls with the number of neurites, and the probability is almost 1 when R is from 3 to 6 (Fig. 8A). Within this range, the probabilities of having a state with 0 or ≥ 2 long neurites is almost zero. At smaller values of R, the probability of having ≥ 2 long neurites increases. In this case the actin waves are distributed among a smaller number of neurites, so that each receives more excitation that can push it across the threshold from short to long. At larger R values the probability of having no long neurites increases, since each neurite receives fewer actin waves and thus it becomes more likely that none will go past the threshold.
When the number of neurites is held constant at R = 4, an optimal range of parameter values exists for either the basal excitation amplitude or the excitation rate (Fig. 8B andC). If either parameter is too small, then the size or frequency of actin waves are too small for any of the neurites to cross over from small to long. If either parameter is too large, then more than one neurite will cross over despite of the negative feedback. The optimum range for all three parameters, R, A 0 , and λ depend on the values of other parameter, as they are determined by the balance among excitatory pulses, retraction, and negative feedback. A change in the value of any one parameter changes the balance.
There have been several modeling studies in which the biophysical mechanism underlying neuronal polarization was competition for a limited supply of some growth factor or structural protein; the neurite acquiring the most becomes an axon [14,16,21,30]. In this section, we show how bistability and excitation amplitude reduction are involved in this mechanism.
To illustrate, we build a simple model based on [16]. Consider a neuron with two neurites whose growth is supported by some growth factor F produced at the cell body. F is transported by actin waves to the neurite tips and diffuses back to the cell body. As in [16], we assume that F slows down the retraction of the neurites. Let C 0 , C 1 and C 2 2.
be the concentrations of F at the cell body and the neurite tips. These quantities, and neurite lengths L i , evolve according to:
The first term on the right-hand side of Eq. ( 13) describes the Fickian flux of F which is proportional to the concentration gradient (C i -C 0 )/L i . Unlike [16], we assume that C 0 is constant, so that the total amount of F is not conserved. (The results below are similar whether the growth factor is conserved or not.) The second term represents actin waves carrying the growth factor F . The term approximates the narrow Guassian spikes used in the original model [16]. Each wave causes a jump of size B in concentration. The first two terms of Eq. ( 14) describe the growth of a neurite that is limited by a common resource that is used up as the neurites become longer (e.g., the protein tubulin, which is a key constituent of microtubules). The third term is the retraction rate that is reduced when growth factor is present. It retains the key properties of the original retraction term in the model in [16]: (1) Concentration at the cell body 2
Table 2: List of parameters in the limited-resource model. We choose a constant excitation amplitude (B = 1), reflecting a constant level of growth factor at the cell body (C 0 = 2). Other parameter values are as in [16].
Using parameter values based on experimental measurements in [16] (see Table 2), Monte Carlo simulations reproduce the emergence of a single axon (Fig. 9A). Although fluctuating, the factor mostly accumulates in the long neurite that ultimately becomes the axon. To see bistability in the model, we replace the pulse term in Eq. ( 13) by its average Bλ/2 to obtain an auxiliary deterministic system:
This system evolves on two disparate time scales, with the growth factor concentrations changing much more rapidly than the neurite lengths. In the quasi-steady state, in which dC i /dt = 0, C i is given by
Substituting this expression for C i into Eq. ( 16), we get a two-dimensional system for L 1 and L 2 , whose phase portrait is shown in Fig. 9B. The two stable equilibria in the phase plane demonstrate the bistability in the model at the two polarized states.
To see that the model employs length-dependent excitation amplitude reduction, we plot K 5 c /(K 5 c + C 5 i ) for the neurite that developed into the axon (i = 1 for the case shown in Fig. 9A), which provides for random fluctuation in its retraction rate. As Fig. 9C shows, the fluctuation becomes progressively smaller as the neurite grows. As a result, the "stochastic noise" in the neurite length is damped as the neurite grows, reflecting an excitation amplitude reduction.
Physically, the decay of the fluctuation results from the decreasing Fickian flux as the neurite grows, which facilitates growth factor accumulation. Length-dependent Fickian flux was also involved in other limited-resource models [21,30].
In the original model [16], the total amount of the growth factor was assumed to be conserved. Therefore, the neurites competed for both the growth factor and structural proteins. Also, the excitation amplitude B was assumed to be proportional to C 0 , which decreased as the neurites grew. This was a second means of reduction in the size of the stochastic noise. Although unnecessary for successful neuronal polarization, as we showed here, these additional mechanisms may help the establishment of a single-axon in a noisy biological environment. Redundancy in biological processes is common in biological systems [24].
In this article, we developed and analyzed a minimal model for achieving neuronal polarization that is based on what we believe to be the two key ingredients of the polarization process: bistability and length-dependent negative feedback. The bistability is necessary for the formation of two distinct classes of neurites (short and long), while the length-dependent negative feedback assures that once a neurite becomes long the others are prohibited from doing so.
While there are several plausible targets of the negative feedback in the minimal model, we demonstrated that one stands out as the most effective in achieving persistent polarization. The success of this mechanism, targeted to the amplitude of stochastic actin waves, was demonstrated in several ways, including the joint probability distribution, Monte Carlo simulations, a large ϵ-committor, a long escape time from a polarized state, and a low escape probability.
Additionally, we found that with this negative feedback mechanism, polarization is more successful if there are more than two neurites competing in the winner-takes-all contest. Finally, we demonstrated that a neuronal polarization model based on competition for a limited growth factor has the same underlying key ingredients as our most successful minimal model: bistability and length-dependent reduction in the excitation magnitude.
The clear distinctions between the axon and other short neurites of a neuron during its polarization indicate an inherent bistability [1,31]. Typically, bistability arises from positive feedback [24], and various sources of positive feedback have been identified. One example involves length-dependent retrograde diffusion flux of polarity effectors [4,16,21]. Another results from the anterograde transportation of polarity effectors that is enhanced by their accumulation at neurite tips [14], possibly due to stabilization of microtubules [22]. Microtubule stabilization was also shown to help the localization of endoplasmic reticulum tubules, which in turn enhanced the stabilization [23]. Some of the signaling pathways involved in polarization are discussed in [10,13,15,20].
Another major element of positive feedback in neurite growth is the autocrine effects of neurotrophic factors such as Brain-Derived Neurotrophic Factor (BDNF) and neurotrophin-3 (NT-3). These factors are released by individual neurites and bind to receptors on the neurites, stimulating their growth. It has been shown that BDNF activation of its receptor TrkB not only promotes neurite growth, but provides positive feedback by promoting BDNF secretion [19]. The impact of the local neurotrophin secretion depends on the receptor density at the neurite tip, and it has been shown that neurotrophin binding to receptors recruits more receptors to the membrane, thus providing positive feedback in the response to the neurotrophin [19,32].
Once an axon has formed, negative feedback mechanisms are necessary to prevent the formation of a second axon. In this study, we examined three different unbiased negative feedback mechanisms. Length-dependent increased retraction prevents the growth of a short neurite by destroying its bistability. This negative feedback may result from long-range signals emitted from neurite tips [33], or from a competition for material proteins [16]. We showed that this mechanism is successful in creating polarized states, but does not prevent flipping between the polarized states, which does not appear to occur in actual neurons. This demonstrates that maintenance of a unique polarized state depends on length-dependent suppression of random actin waves, at least in the case of unbiased negative feedback.
It is certainly possible that some form of biased negative feedback occurs, in which only a long neurite can initiate the negative feedback. One example of this is with the neurotrophin NT-3. This growth factor can accumulate at a long neurite and initiate Ca 2+ waves that travel from the neurite tip back to the cell body, activating the small GTPase RhoA that inhibits growth of all neurites [34]. Thus, a growth factor can contribute to neuronal polarization by both facilitating growth of a neurite exposed to it and by inhibiting the growth of competing neurites.
In addition to increased retraction, we also studied the effects of reducing excitation rate and magnitude. We found that reducing excitation rate was insufficient for preventing the formation of a second axon, while reducing excitation magnitude effectively maintained the polarization of our model neuron. It was observed in previous experiments that the frequency of actin waves (i.e., the average number of actin waves per unit time) and the net growth driven by a single actin wave both decreased after an axon had formed [18]. Our study suggests that the decrease in net growth is a more crucial factor in preventing the formation of a second axon.
Previous studies have proposed a mechanism in which all neurites compete for a limited amount of growth proteins, and the neurite that acquired the most becomes an axon [16,20]. It was assumed that the axon's acquisition of these proteins was facilitated by active anterograde transportation and retrograde diffusion. Using a simplification of one such model [16], we demonstrated that this mechanism exhibits bistability and a length-dependent reduction in excitation magnitude. The excitation magnitude reduction results from the axon's decreased diffusion flux as it grows. Lengthdependent diffusion flux is not the only means of preventing the redistribution of growth proteins. It is also possible that blockage occurs in the long neurites. This was demonstrated in a study showing a novel cytoskeletal mechanism in which a dampened retrograde microtubule network assists in the accumulation of Kinesin-1 in the neurite that becomes the axon [35]. Similar to the effect of slow diffusion, the retrograde transportation of Kinesin-1 is reduced, which prevents its redistribution among all neurites.
Our model's bistability and reduction in excitation magnitude may not completely prevent the emergence of multiple axons, which could be considered a flaw. However, a previous experimental study found that a short neurite was able to develop into an axon when it was mechanically stretched, even after another axon had already formed [31].
Our model easily explains this result, as mechanical stretching can cause a neurite to surpass the threshold length, putting it into the basin of attraction of the higher stable equilibrium, regardless of whether another axon already exists. In contrast, a limited resource model that does not allow for more than one axon could not account for this experimental finding.
In a prior experimental study by Wissner-Gross et al., it was observed that neurons with varying numbers of neurites polarized synchronously [2]. The authors found that prior models based on competition for a limited resource [16,36] failed to replicate this, but instead the polarization time increased with the number of neurites. However, if the amount of the limited resource was increased with the number of neurites, the polarization time was similar for model cells with different numbers of neurites. Indeed, they found that the levels of two polarity factors, Shootin1 and HRas, were both higher in neurons with more neurites. We find similar behavior with our models. In the minimal model (Eq. ( 1)), if the basal excitation amplitude A 0 is properly up-regulated according to the neurite number, the time to polarize will be similar regardless of the number of neurites. Similarly, for our limited-resource model (Eqs. ( 13) and ( 14)), if the concentration of the growth factor at the cell body, C 0 , is adjusted based on the neurite number, the time to polarize will remain unchanged. The Wissner-Gross study also found that the majority of rat hippocampal neurons grown in cell culture had between 5 and 7 neurites [2], suggesting the existence of an optimal range for the number of neurites, as in our Fig. 8, and raising the possibility of a regulatory mechanisms for achieving polarization by modulating both the number of neurites and the levels of effectors.
the previous time step by:
or
To estimate Q (i ′ ,j ′ )→(i,j) , we sample M points uniformly in C (i ′ ,j ′ ) . Starting from each point, we could find a trajectory by solving Eq. ( 1) with a Monte Carlo method. Let M i,j be the number of trajectories that end in C i,j , then M i,j /M approximates Q (i ′ ,j ′ )→(i,j) . This should be repeated many times and the average taken. This is a timeconsuming procedure, so we employ the more efficient procedure developed in [45]. Starting from each sample point (i ′ , j ′ ) we solve Eq. ( 1) without the stochastic term over time ∆t:
i,j be the number of trajectories that end in C i,j , then M
i,j /M approximates the probability of transition from C (i ′ ,j ′ ) to C i,j , provided no pulse occurs during ∆t. (The superscript "(d)" represents "deterministic".) Then we consider the case where a single pulse occurs during ∆t. Since we describe actin waves as a Poisson process, the time of the occurrence is uniformly distributed within ∆t [46]. Thus, we solve Eq. ( 23) over ∆t/2, then randomly choose a length from L 1 and L 2 and add A to it, and finally solving Eq. ( 23) over the rest of the time interval ∆t/2 (Fig. 10).
The result is a trajectory with a single discontinuity. We repeat the same calculation for all sample points. Let M (s) i,j be the number of trajectories that end in C i,j , then M
i,j /M approximates the probability of transition from C (i ′ ,j ′ )
to C i,j , provided a single pulse occurs during this short period. (The superscript "(s)" represents "stochastic".) Since the probability of having two or more pulses during ∆t is of O(∆t 2 ), we neglect this probability and approximate Q (i ′ ,j ′ )→(i,j) as
where λ∆t is the first order approximation of the probability of having a single pulse during ∆t. In principle, one could refine the approximation by dividing ∆t into more subintervals.
Given an initial distribution p(0), we can calculate p(n) iteratively with Eq. (21). To obtain p(0), suppose that the system starts from (L 1 (0), L 2 (0)). We find the cell C i,j containing this point and set the corresponding p i,j (0) to be 1 and all other probabilities to be 0. To estimate the limiting distribution p(∞), we iterate according to Eq. ( 21), until the change in p(n) becomes negligible. To speed up the iteration, we utilize Q 2 k = (Q 2 (k-1) ) 2 , such that p(2 k ) can be obtained with k iterations.
In addition to solving the distribution, we will also use the GCM to calculate various first passage probabilities and mean first passage times. This requires modification of the transition matrix Q. Suppose we are interested in finding the probability that the two-neurite system enters a specific region D and the mean entering time. For a cell centered within D, the transition probability Q (i ′ ,j ′ )→(i,j) is modified as Q(i ′ ,j ′ )→(i,j) :
This makes D an absorbing region, which means that once the system enters D, it is frozen and cannot make further transitions. Let the modified transition matrix be Q and the resulting distribution be p(n). The probability of entering formed by all q i,j 's is called the ϵ-committor, namely q = (q 1,1 , q 1,2 , . . . , q N,N ) T .
Since the probability of being absorbed at each time step by either Z 1 or Z 2 is ϵ, the mean time till absorption is ∆t/ϵ, where ∆t is the step size used in the GCM. Over such a timescale, if the system starts from C i,j and spends most of the time within R S , it will have a high probability of being absorbed into Z 1 , i.e., q i,j will be close to 1. Conversely, if the system never enters R S or quickly leaves it without coming back, q i,j will be close to 0 [29]. Therefore, q i,j characterizes the attracting strength of the region R S over a timescale of ∆t/ϵ, when the system starts from C i,j . By choosing a starting cell close to the top left peak of a bimodal distribution and changing the value of ϵ, the resulting q i,j quantifies the persistence of the polarized state S over different timescales under the corresponding negative feedback mechanism. Specifically, we choose the cell at [0, 6] when the increased retraction is implemented, and the cell at [0,7] when the excitation rate or amplitude reduction is implemented. The corresponding probability q i,j (R S , ϵ) is denoted by C ϵ for notational simplicity.
To calculate q use the formula in [29], given as
where I is an N 2 -by-N 2 identity matrix and Q the transition matrix given by Eq. (20). I R S is an N 2 -by-1 indicator vector defined as
where