RETRACTED: A mathematical model explains saturating axon guidance responses to molecular gradients
Abstract
Correct wiring is crucial for the proper functioning of the nervous system. Molecular gradients provide critical signals to guide growth cones, which are the motile tips of developing axons, to their targets. However, in vitro, growth cones trace highly stochastic trajectories, and exactly how molecular gradients bias their movement is unclear. Here, we introduce a mathematical model based on persistence, bias, and noise to describe this behaviour, constrained directly by measurements of the detailed statistics of growth cone movements in both attractive and repulsive gradients in a microfluidic device. This model provides a mathematical explanation for why average axon turning angles in gradients in vitro saturate very rapidly with time at relatively small values. This work introduces the most accurate predictive model of growth cone trajectories to date, and deepens our understanding of axon guidance events both in vitro and in vivo.
https://doi.org/10.7554/eLife.12248.001eLife digest
For your brain to work, millions of nerve cells have to be connected together precisely. To achieve this, growing nerve fibres navigate through the developing brain by following chemical cues. One important such cue is how the concentration of a chemical varies with distance across the brain. This variation is known as a chemical gradient. However we still don't fully understand exactly how nerve fibres use such gradients to steer themselves.
Nguyen et al. have built a mathematical model that accounts for the way that nerve fibres respond to chemical gradients, and showed that the model closely matched new experimental data on the growth of nerve fibres from the rat brain. The model implies that, under some conditions, nerve fibres turn surprisingly little in response to a chemical gradient.
Nguyen et al.’s model can now be used to predict nerve fibre responses in many other situations, and could help researchers to understand more about how the brain becomes wired up during development. The model could also reveal more about the conditions that are needed to cause nerve fibres to turn sharply in response to chemical gradients.
https://doi.org/10.7554/eLife.12248.002Introduction
For the brain to function correctly, it must be wired correctly. Indeed, many neurodevelopmental disorders are likely the result of wiring defects (Yaron and Zheng, 2007; Geschwind and Levitt, 2007; Lin et al., 2009; Stoeckli, 2012). Axon guidance, where axons grow and navigate to their targets, occurs primarily via the sensing of molecular cues in the environment. A critical mechanism by which such cues act is believed to be concentration gradients, causing axons to be attracted or repelled in particular directions (Mortimer et al., 2008; Lowery and Van Vactor, 2009). However, despite major advances in understanding which molecules are involved in this process (TessierLavigne and Goodman, 1996; Dickson, 2002; Chilton, 2006; Kolodkin and TessierLavigne, 2011), an accurate quantitative model describing how axon trajectories are influenced by such guidance cues is still lacking.
In vivo, axon trajectories may potentially be influenced by many cues. In vitro assays allow individual influences, such as that from the concentration gradient of a single guidance factor, to be isolated and quantified. A substantial mystery posed by in vitro axonal chemotaxis assays is the relatively weak turning produced, even over long periods of time. The naive prediction that axons would promptly turn until they become fully aligned with the gradient turns out not to be true. In an early study of chemotactic responses of chick sensory neurons to a gradient of nerve growth factor in a diffusion chamber, only 60% of nerve tips were preferentially directed toward the gradient direction after 46 hr of growth (Letourneau, 1978). The growth cone turning assay over 1–2 hr produces average turning angles typically ranging from 10 to 25°, with high variability (Song et al., 1997; Höpker et al., 1999; Xiang et al., 2002; Ming et al., 2002; Thompson et al., 2011). A similarly weak response is observed in the Dunn chamber (Kent et al., 2010; Dudanova et al., 2012; Ruiz de Almodovar et al., 2011; Dudanova et al., 2010; Yam et al., 2009). More recent studies using microfluidic technologies over timescales ranging from hours to days have also elicited average axon turning angles only up to 10–15° (Wang et al., 2008; Morel et al., 2012; Taylor et al., 2015; Sloan et al., 2015). Why average turning angles are so small, and what this means for axon guidance in vivo, are unclear.
One of the key properties of in vitro axon growth that might explain this mystery is that it is often very straight (Katz et al., 1984; Katz, 1985). Axons are under mechanical tension from the pull of the growth cone (Bray, 1973; Bray, 1979), and this tension stimulates the elongation of the axon by stretching (Bray, 1984; Zheng et al., 1991). Traction forces generated in the growth cone arise from the coupling of the continuous retrograde flow of actin to the substrate through adhesion receptors (Franze and Guck, 2010; Betz et al., 2011; Athamneh and Suter, 2015). For reasons which are not clear, axons tend not to bend and follow the highly random movements of their growth cones. Rather, they usually form a straight line between their tip and a location where they are firmly attached to the substrate (i.e. a focal adhesion (Kaverina et al., 2002)). We call such locations anchor points; they can be at the soma, at a branch point, or at some other seemingly sporadic location along the axon. Although it is not clear how this tension leads to elongation, the growth cone advances largely in the stretch direction along the axon, resulting in relatively straight paths.
To determine quantitatively what effect this might have on axonal trajectories requires mathematical modelling. Growth cone movements were first analyzed in detail in (Katz et al., 1984; Katz, 1985). Subsequently, various phenomenological models have been built that differ as to how they treat stochasticity, and mechanisms for directional preference, namely turning or growth rate modulation. One set attempted to fit the dynamics of growth cone movement to a random walk with drift (Buettner et al., 1994; Odde and Buettner, 1995; Maskery et al., 2004; Pearson et al., 2011). Li et al. simulated trajectories by assuming the turning angle of the growth cone is in proportion to the angle between the neurite and the resultant filopodial tension (Li et al., 1995). In (Borisyuk et al., 2008), the axon growth angle depends on the tendency to turn toward the gradient angle and noise. The noise term is small (2–5°), leading to straight paths that resemble axon growth in the tadpole spinal cord. Another set of models has concentrated on how asymmetric receptor binding across the growth cone might be used as the basis of a turning signal (Goodhill et al., 2004; Xu et al., 2005; Mortimer et al., 2010), but without considering the consequence for whole trajectories. A third group of models considers the possibility that the velocity of the growth cone is influenced by an attractive gradient from the target cells, and chemoattractants and chemorepellants released from other growth cones and itself (Hentschel and Ooyen, 1999; Krottje and van Ooyen, 2007). However, none of these models has been closely compared with the details of experimentally measured trajectories in gradients, and parameters such as variability in step sizes, the distribution of instantaneous turning angles, and straightness of real paths, have not been addressed. Thus, the question of whether there is a model that can adequately capture all these characteristics of real trajectories remains open. Without such a model, it is difficult to determine if trajectories observed in vivo are in fact consistent with gradient guidance.
Here, we present a new computational model for axonal trajectories based on the combined influence of anchor points, a tendency to turn toward the gradient direction, and random noise. We found experimentally that the gradient had no effect on the step sizes; thus, we only model the turning angles. Critically, the model predicts rapid near saturation of average turning angles with time. To test this model quantitatively, we then introduce a new microfluidics assay for studying axonal response to gradients, and using timelapse imaging characterize the behavior of axons over several hours of growth in both attractive and repulsive gradients. We find that our model fits the behavior observed very closely. We then investigate by simulation the effect of increasing the number of anchor points, and find that this increases the average fidelity of turning but at the cost of higher variability. Together, this work both explains why turning response to gradient saturates so rapidly and reveals the quantitative principles that are required to reproduce accurately in vitro axonal trajectories in response to chemotactic gradients. The model identifies straightness as a limiting factor on how much axons can turn and suggests that the frequency of anchor points plays a key role in the axonal turning response to a gradient.
Results
A mathematical model of growth cone trajectories
We modelled three basic influences on the direction of axon growth: a tendency to grow straight, the effect of a chemotactic gradient, if present, and random movement noise. In a fixed coordinate system with arbitrary zero angle direction, we define θ(t) as the bearing of the growth cone at time step $t$, ϕ($t$) as the angle of the vector connecting the growth cone to its anchor point, and Ψ as the gradient direction (terminology is summarized in Table 1). We define ‘bearing change’ as Δθ($t$), the change in θ(t) at time step $t$, distinct from ‘turning angle’ ψ_{turn}, the total change in θ from the initial direction of growth over long periods of time. For simulations we identify each timestep as 5 min of real time. The model (Figure 1) is then

Figure 1—source code 1
 https://doi.org/10.7554/eLife.12248.005
where $a$ scales persistence to move in the same direction as the overall direction of the axon, $b$ scales the bias due to the gradient, and ξ is random noise in the bearing changes. The symbol ∠($x,y$) denotes the signed angle between the unit vectors with angles $x$ and $y$, and constrains the resultant angle to be between −π and π. The step size is the distance moved after one time step and will later be estimated empirically.
We consider first the noiseless case (ξ = 0) in long and finitetime regimes, and then consider the effects of noise. Figure 1C shows the results of setting ξ = 0, with a fixed step size of $s$= 3 μm, and simulating the model for long times with the same $a$ = 1 and different values of $b$ (0.1, 0.2 and 0.3). Turning angles rapidly saturate, which can be understood analytically (see 'Materials and methods'): in the $t\to \infty $ limit, the growth cone angle follows a power law with respect to time $\varphi \left(t\right)\propto {t}^{(\alpha 1)}$ or $\mathrm{log}\left(\varphi \right(t\left)\right)=const+(\alpha 1)\mathrm{log}\left(t\right)$ (where $\alpha =a/(a+b)$) (Figure 1C). This relationship generally holds for $t>\mathrm{exp}\left(4\right)\approx 4$ h, meaning that for long times, the rate of change of angle decreases and this rate is determined by the power law exponent $b$/($a$ + $b$). Since comparison with empirical data (see later) shows that the biologically relevant regime is $b$ ≪ $a$, the exponent is generally small. Thus, while ultimately axons in the model do eventually align with the gradient, this process takes an exceedingly long time. This explains the slow and decreasing change in the turning angle over time in the noiseless case.
The finite $t$ regime of this equation is difficult to solve analytically, since $\varphi \left(t\right)$ depends on the entire history of growth cone movements. Simulations using different combinations of $a$ and $b$ are shown in Figure 1D. For the cases of $a\ne 0$, after 150 time steps (12.5 hr of real time), the resultant turning angle was far from completely aligned with the gradient. Although the bias term bent the trajectory in the direction of the gradient, there was a straightening effect due to the persistence term, constantly pulling the growth cone toward the overall growth direction of the axon. As expected, the pull due to the gradient increased with larger b (Figure 1D). Thus, the persistence term prevented the axon from completely aligning with the gradient. Also apparent is that without noise, the trajectories were all very straight (with straightness index [see 'Materials and methods'] greater than 0.98). Thus, the microscopic constraint imposed by the persistence term leads to the macroscopic phenomenon of incomplete turning.
When we introduced Gaussian noise into the bearing changes (in radians) $(\xi ~N(0,\sigma ),\sigma =\pi /4)$ with the same parameters and initial conditions as above for 1000 axons, the behavior was qualitatively similar: after an initial period of relatively rapid turning, turning angles tended to an almost steady state which was not aligned with the gradient even after a long time. However, the average final turning angle was even less than that of the noiseless case. This is because the noise created more random wandering of the growth cone, further reducing the directional effect of the gradient (Figure 2A–B). After 20–40 min, for $b$ = 0.1, the turning angle distribution of the population was 7 ± 25° (mean ± std). Assuming a normal distribution of turning angles, this means that many of the axons were no longer roughly perpendicular to the gradient, and thus only continued to turn extremely slowly. Therefore, over time the influence of the gradient on the whole population of axons weakened. The persistence term also created a resistance against large turns due to the gradient. Increasing $b/a$ increased the turning angle, but did not alter its rapid saturation with time (Figure 2C). Lastly, we examined the effect on the straightness by varying the standand deviation of the noise σ from 0 to π radians. As the steps became more noisy, the paths became less straight (Figure 2D).

Figure 2—source code 1
 https://doi.org/10.7554/eLife.12248.007
In summary, the noiseless case generated very straight axons and growth cone angles that followed power laws with respect to time in the long time limit. Similarly, in the noisy case, the rate of change of the average turning angle was initially rapid and then slowed down even more rapidly with time. In both cases, the persistence term was a limiting factor on how much and how fast the axons could turn. Thus, this model captures, at least qualitatively, the behavior that axons turn only slightly in gradients, and even for long times do not generally become completely aligned with the gradient.
Stable gradient generation for guidance assays
Having established the basic behaviour of the model, we then asked whether it could reproduce in detail real axon trajectory statistics. We therefore analyzed the trajectories of superior cervical ganglion neurons in a new microfluidics device (Figure 3A–C). This device generated linear gradients, by the mixing of high and low solutions of a chemotropic factor. The gradient was visualized using 40 kDa dextrantetramethylrhodamine (Figure 3D), and gradients were stable for at least 20 hr (Figure 3E,F).

Figure 3—source data 1
 https://doi.org/10.7554/eLife.12248.009
SCG neurons were guided in the microfluidic assay
We measured the response to nerve growth factor (NGF) gradients of axons from dissociated P1P3 SCG neurons. We chose this model system because almost 100% of these neurons express the NGF receptor TrkA (Wetmore and Olson, 1995; Verge et al., 1992).
Three conditions were investigated: a control without flow or gradient, an attractive gradient of nerve growth factor (NGF), and a gradient of NGF with added KT5720, which converts attraction to repulsion by lowering levels of cAMP in the growth cone (Song et al., 1997). Cells were injected into the growing chambers and grown for 2 hr before gradient onset. In the control condition, cells were grown over several hours with 0.3 nM NGF. In the NGF gradient condition, two solutions of concentrations 0 nM and 10 nM NGF were pumped into the growing chamber through the two inlets. Previous work using Scatchard analysis estimated that K_{d} = 0.9 ± 0.3 nM (Wehrman et al., 2007) and showed that SCG neuronal outgrowth is severely inhibited at the saturating NGF concentration of 40 nM (Ohta et al., 1990). Given the healthy growth in our assay, it is clear that the concentration in the gradient condition was below saturation point. We analyzed trajectories for 300 axons per condition. These were obtained from 23 individual chambers in the control case, 27 chambers in the NGF gradient case, and 24 chambers in the NGF gradient plus KT5720 case. In most experiments, 2 chambers were run in parallel, so the total numbers of experiments in each case were 12, 15, and 13, respectively. For the measurement of turning angles, we selected only axons that started growing between 70° and 110° relative to the gradient (when present). An asymmetric concentration field of guidance cue across the growth cone leads to turning (Song et al., 1997; Hentschel and Ooyen, 1999; Xiang et al., 2002; Ming et al., 2002) and axons growing in this range experienced between 94% (i.e. sin 70°) to 100% (i.e. sin 90°) of the maximum possible concentration difference across the growth cone. Thus, we expected the impact of the gradient would be strongest on these axons (Figure 4A). We tracked the growth cones every 5 min for as long as possible until they collapsed or branched or collided with other cells, axons or the edges of the chamber. The SCG axons were clearly attracted in the NGF gradient (Figure 4B). When the protein kinase A inhibitor KT5720 was added to the high and low solutions at concentration 70 nM, attraction was converted into repulsion as previously described (Song et al., 1997; Forbes et al., 2012) (Figure 4B). These results confirm that the gradient in the microfluidic assay elicited a guidance response in SCG axons. From the timelapse imaging data, we then selected the subset of axons that did not branch or retract following growth for several hours in the attractive NGF gradient, and measured the turning angles of the population over time. The average turning angle reached the steady state quickly and did not increase significantly with time, matching the prediction of the model (Figure 4C).

Figure 4—source data 1
 https://doi.org/10.7554/eLife.12248.011
The gradient did not affect axon branching
One possible way that the gradient could affect the axons is by causing biased branching [c.f. (Simpson et al., 2013)], or changes in branching rates. To test whether the NGF gradient changed the branch extension and retraction rates, we compared the number of branches per cell after 5 hr of growth and did not detect any difference (p = 0.9 Kolmogorov–Smirnov test, Figure 5A). We measured the intervals between successive branching events in the same cell in each condition and did not find any difference in the branching rate (p = 0.7 Kolmogorov–Smirnov test, Figure 5B). Similarly, the lifetimes of the branches were unaffected by the gradient (p = 0.2 Kolmogorov–Smirnov test, Figure 5C). We counted the number of branches pointing up and down the gradient per cell and did not find any difference (p = 0.8 Kolmogorov–Smirnov test, Figure 5D). Thus, the gradient had no effect on axon branching and retraction.

Figure 5—source data 1
 https://doi.org/10.7554/eLife.12248.013
Flow did not affect the statistics of steps
To test whether fluid flow in the chamber biased the statistics of the steps, axons growing in the gradient condition with fluid flow were divided into four quadrants with different relative angles to the fluid flow: two quadrants growing perpendicular to the flow, one quadrant growing with the flow, and the other growing against the flow (Figure 6A). Comparing the distribution of bearing changes between the four quandrants, and with axons from the control condition without any flow, showed no influence of the flow (p = 0.7 in Figure 6B and p = 0.4 in Figure 6C, KruskalWallis test). The means of the bearing changes in quadrants 2 and 4 were nonzero, and the bearing changes accumulated over time to result in a nonzero average turning angle of the population. However, these differences in means were very small (approximately 1°), and there were no significant statistical difference among the distributions. Therefore, the positive turning angles in the NGF gradient (and negative turning angles in the NGF gradient + KT5720) were due to the effect of the gradient, not bias from the flow.

Figure 6—source data 1
 https://doi.org/10.7554/eLife.12248.015
Growth cone trajectories were generally straight
Axon growth is shown in Videos 1–3 and Figure 7A–C. The growth cone often wandered quite randomly, but nevertheless usually the axon segment remained very straight from the growth cone to the cell body or last axon branch point (Figure 7A–C). This implies that often the entire axon segment was pulled sideways across the substrate (as can be seen directly in the movies). Thus, despite the irregular trajectory of the centre of the growth cone, the tension force on the growth cone from the axon was usually pointing directly back to the last anchor point, consistent with the assumptions of the model. To quantify this further, we measured the angle of this 20 μm segment and the angle to the anchor point ϕ and found them to be almost the same (Figure 7D,E). Thus, we can understand the term ϕ as the tension due to the most distal segment.

Figure 7—source data 1
 https://doi.org/10.7554/eLife.12248.017
The trajectories (i.e. the locus of the centre of the growth cone) in three conditions are plotted in Figures 8–10. Note that these paths are not the same as the final image of the axon, which generally pointed straight back from the final position of the growth cone to the anchor point. Visually, the paths appear mostly straight with occasional large turns, consistent with a long tail for the bearing change distribution. The mean straightness index for the trajectories was S = 0.72 (Figure 11A).

Figure 11—source data 1
 https://doi.org/10.7554/eLife.12248.025
Step size and bearing change distributions were similar across conditions
There was little correlation between the bearing change magnitude and step size (Figure 11B). The distribution of bearing changes in radians was well fitted by a mixture of a von Mises and a uniform distribution (−π < x < π) (Figure 11C). That is there was a great deal of randomness in bearing changes, but with a peak in probability near the forward direction. Thus, growth cones tended to move in a straight line instead of turning uniformly randomly. This is inconsistent with the assumptions of several previous models (Buettner et al., 1994; Odde and Buettner, 1995; Maskery et al., 2004).
Accumulating across all the growth cones, the distributions of step sizes over 5 min were statistically indistinguishable across the three conditions (KruskalWallis test p = 0.35), and were well fitted by a gamma distribution (Figure 11D). That is, the most likely step size was around 0.5 μm/min, but the distribution had a long tail, so that longer step sizes were also seen. The distribution of step sizes for each individual growth cone were also well fit by gamma distributions (Figure 11E). However, individual growth cones had idiosyncratic mean values. The distribution of these mean values could be well fitted by a Gaussian distribution (Figure 11F).
Nevertheless, the mean square displacement was clearly not linear, implying that a simple random walk is not suitable to describe the movement (Figure 11G). Successive steps were anticorrelated (Figure 11H), which was not accounted for in a previous model (Pearson et al., 2011). This helps the paths remain relatively straight: if successive steps were positively correlated, the paths would become more bent over time. Due to large noise in the bearing changes, bearing changes more than one step apart were uncorrelated.
Turning angles over time were well predicted by the model
Having established the key statistics of steps from the data, we now asked if the simple model in Equation (1) could replicate the observed trajectories and explain the phenomenon of saturated turning. We sampled the mean speed ${v}_{mean}$ of each growth cone from a truncated Gaussian distribution of mean 0.7 μm/min and standard deviation 0.24 μm/min. At each time point (5min interval), the growth cone sampled a step size from the gamma distribution $\Gamma \left(4\u2215u,{v}_{mean}*u\u22154\right)$ where u was a uniform random number. The bearing changes evolved according to Equation (1). We found that the random noise ξ in bearing changes (in radians) could be well fit by the mixture von Mises distribution
where $c$ and $d$ are parameters to be fit. This distribution is not necessarily the same as that of the bearing changes in Figure 11F. As the bearing change is the sum of three random terms, its distribution is broader than the distribution of the noise term. To estimate the four free parameters $a,b,c,d$, we input the initiation angles $\varphi \left(0\right)$ and used the model to generate the distribution of turning angles ${\psi}_{turn}$. We then estimated the likelihood function that the turning angle data was generated from the model with the given parameters $L\left({\psi}_{turn}\righta,b,c,d,\varphi \left(0\right))$ .
We found the values of $a,b,c,d$ that maximized the likelihood (the parameters mostly likely to have generated our empirical turning angle data) were $a$ = 0.7, $b$ = 0.09, $c$ = 0.75, $d$ = 6. The statistics over 5000 simulated trajectories using these parameters are shown in Figure 12, and some example trajectories are shown in Figure 13. Similarly, we fitted the model to the control data and repulsive gradient and found $b$ = 0.002 and $b$=−0.08, respectively, with other parameters remaining at very similar values as before. Notably, the simulated turning angles changed little over time (Figure 12A), consistent with our preliminary prediction and the data (Figures 1E, 4C). Thus, with realistic step sizes and bearing change noise distributions, the model was able to capture the phenomenon of saturated turning with quantitatively accurate means and variances over time. The distribution of turning angles and step sizes also closely matched the real data (Figure 12B,C). The simulated straightness distribution was also very smilar to the real distribution (p = 0.2, ttest) (Figure 12D). The model also captured the distribution of bearing changes, the mean square displacement and the anticorrelation between successive bearing change, which was a consequence of the persistence term straightening the paths (Figure 12D–F). If successive steps were positively correlated, the paths would become more bent over time. This correlation was rapidly lost beyond one time lag because of the large noise.

Figure 12—source code 1
 https://doi.org/10.7554/eLife.12248.027
Unlike previous models, we did not assume constant steps or a uniform distribution of bearing changes but rather derived these from empirical data. The model was then able to predict the evolution of the average turning angle over time, the straightness profile and the anticorrelation in bearing changes. Most importantly, it could explain the phenomenon of slow and saturated turning, due to a weak bias term relative to the persistence term. A microscopic factor in each step led to a macroscopic phenomenon of limited, variable turning and straight paths. This often overlooked feature of axon growth turned out to be critical in our model in limiting the overall turning. We also found little difference between the attractive and repulsive case, indicating that attractive and repulsive gradients employed similar mechanisms and could not reduce the variability of axon trajectories.
Multiple anchor points achieved sharp turns but also increased variability
The in vitro data we have presented here was wellfitted by assuming the only anchor point is where the axon emerges from the soma or the branch point. However, the in vivo environment is much more complex, and axons may establish anchor points with the substrate at multiple positions as they extend. We therefore investigated in the model what effect this would have on turning angles. We assumed that at each timestep, the probability of that point becoming a new anchor point was fixed, while leaving the evolution at each step as before. The average number of anchor points per timestep (i.e. 5 min) is denoted by r. We analyzed two cases: anchoring probabilistically at each time step (Figure 14A–C), and anchoring at regular intervals (Figure 14D–F). We simulated the trajectories for T = 150 timesteps with the same parameters as Figure 2A ($a$ = 1, $b$ = 0.1, ξ = 0 or ξ ∼ N (0, π/4) radians). In both cases, more anchor points led to sharp turns in the trajectories and larger mean turning angle (Figure 14G), since the growth cone now was more free from its initial position. However, it increased the variability in the turning (Figure 14G). Given the same rate of anchoring, whether the growth cone put down new anchor points probabilistically or regularly made little difference to the mean turning. We compared the mean square final angle $\u27e8\varphi {\left(T\right)}^{2}\u27e9$ which is the sum of the bias ${\u27e8\varphi \left(T\right)\u27e9}^{2}$ and the variance $\left(\varphi \left(T\right)\right)$. Ideally, the growth cone should completely align with the gradient, that is $\varphi \left(T\right)=0$. Figure 14H shows the bias/variance trade off. Although more anchor points introduced larger variance in the final angle, they achieved greater turning, that is smaller $\varphi \left(T\right)$. Figure 14I shows that in the case without movement noise (ξ = 0) and regular anchoring, the growth cone angle $\varphi \left(t\right)$ also followed a power law in the large t limit, similarly to the case without anchoring (Figure 1C). However, the exponent of this power law was larger, demonstrated by a steeper slope between log($t$) and log(ϕ($t$)), meaning that the growth cone aligned with the gradient faster with more anchor points. Thus, increasing the rate of anchoring leads to stronger turning, but increases the variance of the responses.

Figure 14—source code 1
 https://doi.org/10.7554/eLife.12248.030
Discussion
Here, we presented a model of axon trajectories in gradients and helped resolve the mystery of why axon turning angles in gradients saturate over time in vitro, revealing an important factor limiting axon turning. We found that the movement of the growth cone was strongly influenced by the axon’s tendency to maintain a straight trajectory forward, limiting the directional effect of the gradient and preventing the axon from aligning with the gradient even after a long time. Our model predicted that, averaged over a large population of axons, the initial rate of turning drops rapidly over a short period of time (20–40 min). The model shows that adding more anchor points can give the growth cone more flexibility and produce larger average turning, but also increases the variability. Thus, we predict that different substrates, producing different densities of anchor points, could result in different trajectories for the same the gradient conditions.
The application of forces to axons can induce rapid elongation without axonal thinning, and thus stretch can stimulate growth (Suter and Miller, 2011). Furthermore, stretch can also regulate the mode of growth. When axons are tightly bound to a sticky substrate, stretching only happens at the tip and axons elongate by tip growth. In contrast, if axons grow relatively unattached to the substrate, they will lengthen by stretching due to the pull of the growth cone (Chang et al., 1998; O'Toole et al., 2008), which appears to be the case in our experimental condition. The tension along the axon will cause stimulation of growth in the existing direction producing straighter trajectories. The stiffness of axons is also important (Rajagopalan et al., 2010), and stiffer axons will likely have higher persistence due to their more limited ability to bend.
This tension results from cytoskeletal coupling with adhesive interactions to the substrate and is critical to growth cone migration (Heidemann et al., 1997; Spire, 2009). Although anchor points are an abstraction in our model, their biological implementation may be focal adhesions. Only at these points is the axon firmly attached to the substrate. There is a number of ways in which anchor points could be investigated experimentally in future work. Axons could be stained for proteins such as integrins (Kaverina et al., 2002) to test whether their distribution is strongly localized to particular points along the axon. We also predict that applying force orthogonal to the direction of axon growth, for instance by using a pipette to puff liquid at different locations along the axon, would cause a deflection of the axon of a size related to the distance from the nearest anchor point. A similar experiment was performed using a glass needle to tow axons (O'Toole et al., 2008). It was observed that the distal region of the axon was free of the substrate, while the proximal region was firmly attached. In addition, it could be possible to determine the internal stress field of an axon, as has been done for growth cones (Betz et al., 2011): we would expect the stress to in general be different on the two sides of an anchor point. The density of anchor points will depend on the components of the extracellular matrix (ECM). In our experiments, on laminin, they appeared to be rare. This might be because adhesion points are expensive to produce and the axon can grow faster when it is not attached to the substrate (Chang et al., 1998). However, the biological factors governing when new anchor points are generated are unknown.
Tension is also dependent on cell type and two main properties of the substrate: stiffness and ECM components. Our data comes from peripheral nervous system (PNS) neurons growing on a laminin substrate that is hard rather than gellike, and other cell types of different substrates might have different behaviours. Central nervous system (CNS) and PNS neurons have different sensitivities to substrate stiffness due to adaptation to their natural environments (Koch et al., 2012), and traction force in vitro increased on stiffer substrates (Koch et al., 2012). Substrates with different ECM components differentially promote growth cone motility and point contact formation. For example, growth cones are more highly motile and neurites extend more rapidly on laminin than fibronectin because point contacts have higher turnover rate (Robles and Gomez, 2006).
Overall, our work suggests that without many anchor points, cues additional to gradients may be necessary for axons to reliably find their targets in vivo (unless the motility noise is for some reason much lower in vivo than in vitro). These could include mechanical cues and axonaxon interactions. To understand such interactions, it is important to generate assays with realistic substrates suitable for different cell types. Recent 3D culture models, in which cells are grown with a protein scaffold, can capture some aspects of the tissue environments instead of hard surfaces (Cullen et al., 2007). It will be interesting to see how different ECM properties lead to changes in trajectories and whether they can facilitate more reliable turning.
In conclusion, we have presented a simple mathematical model which gives accurate quantitative predictions of the properties of axonal trajectories in a microfluidicsbased in vitro gradient assay. The model identifies the key importance of anchor points in controlling turning and provides an explanation for why axonal turning in gradients in vitro tend to saturate rapidly at small turning angles. This model provides a predictive framework which can be used to test whether axonal trajectories observed in vivo can be explained purely in terms of gradient guidance, or whether additional guidance mechanisms are also required.
Materials and methods
Microfluidics chamber fabrication
Request a detailed protocolThe device was designed in AutoCAD 2013 based on the pattern in (Dertinger et al., 2001). The microfluidic master molds were prepared by standard soft photolithography techniques. The silicon wafer (M.M.R.C. Pty Ltd) was coated with a 50 μm thick layer of SU8 2100 relief (MicroChem, Westborough, MA) with a spincoater (EVG). The master was printed onto highprecision photoplates (KonicaMinolta) with a photoplotter (EVG). Before replica moulding, the silicon master was silanized with vapour phase silane, under vacuum for 1 hr at room temperature to prevent polydimethylsiloxane (PDMS) adhering to the master.
PDMS base elastomer (Sylgard 184, Dow Corning, Midland, MI) and silicon elastomer curing agent in a 10:1 (m/m) ratio were thoroughly mixed and poured on the master to a depth of about 4 mm. The mold with PDMS was degassed in a vacuum chamber for 2 hr and baked at 80°C for 2 hr. The PDMS replica was then peeled from the silicon wafer and cut into individual chambers. Holes were cored in the PDMS chambers using a 0.75 mm corer (Harris UniCore). To bond the PDMS chamber to a plastic tissueculture petri dish, the dish and the chamber were plasma treated at 100W at a pressure of 380 mTorr for 40 s. The plate was then covered with (3Aminopropyl)triethoxysilane (APTES, SigmaAldrich) solution (5% APTES in 70% ethanol). The plastic dish was washed thoroughly with water, air dried, and the chamber was pressed onto the dish. To avoid air bubbles forming, the dish was filled with distilled water and degassed in the vacuum chamber for 10 min before use.
Primary superior cervical ganglion (SCG) cell culture
Request a detailed protocolSCG neurons were harvested from postnatal day 1–3 Wistar rat pups. The SCGs were then cut into thirds, incubated in 0.25% trypsin (GIBCO) at 37°C for 30 min and then triturated through a flamedpolished Pasteur pipette for 5 min to dissociate individual cells. The growth solution was OptiMEM solution (GIBCO) containing 1X penicillin/streptomycin, 10 μg/mL mouse laminin, 4% (v/v) fetal calf serum, 2% B27 supplement (Life Technologies), and 0.3 nM NGF. The cells were suspended in the growth solution and injected in the microfluidic chamber using a 100μL glass syringe (SGE Analytical Science).
Gradient generation and measurement
Request a detailed protocolTwo syringes were filled with either the high (10 nM) or low (0 nM) NGF solution. The high solution contained 0.1% (v/v) 40 kDadextran fluorescently labelled with tetramethylrhodamine to visualize the gradient. After the cells were seeded, the syringes containing the high and low solutions were connected to the chamber using polyethylene tubings. The chamber was moved to an incubated inverted microscope (Zeiss AxioObserver). The syringes were attached to a Harvard PHD pump and the flow rate was set at 10 μL/hr. After the flow had been started, fluorescent images of the chamber were taken in Zen software. Background intensity outside the chamber was subtracted from the images. The average brightness intensity and variations over time were then calculated across the chamber. A gradient of fluorescence confirmed that the NGF gradient had formed in the growth chamber. To generate a repulsive gradient, KT5720 (Alexis Biochemicals), a specific inhibitor of protein kinase A (PKA), was added into both the high and low solutions at a concentration of 70 nM.
Tracking growth cone trajectories
Request a detailed protocolAfter the onset of the gradient, the axons were imaged every 5 min for 6 hr using Zeiss Zen software. After data acquisition, axons of 30 μm length, growing in all directions, that did not branch or retract in at least 80 min, were chosen for measurements. All axons were tracked manually using customized MATLAB software (The MathWorks) for as long as possible until they branched or retracted. A 5min time interval was chosen because, for smaller intervals, variability in identifying the centre of the growth cone was larger than the net movement between frames. The point where the axon attaches to the cell body or the main branch was considered the anchor point.
The straightness index
Request a detailed protocolThe straightness index S is the inverse of tortuosity, and compares the overall net displacement G of a path with the total path length T (Codling et al., 2008). Consider a walk that starts at location $\left({x}_{0},{y}_{0}\right)$, and after n steps of lengths l_{j} (j = 1...n) finishes at $\left({x}_{n},{y}_{n}\right)$. The straightness index is given by:
This index is between 0 and 1, where 1 corresponds to movement in a straight line and 0 corresponds to a walk that returns to the origin. The closer this index is to 1, the straighter the trajectory is. Obviously, $S$ depends on the time interval used for tracing but can be used to compare conditions, which all have the same time interval.
Modeling growth cone trajectories
Request a detailed protocolAll parameters of the model are summarized in Table 1. We consider a model which is a discretized random walk in which we separate the length and directions of the steps (Figure 1A). We discretized the axons at a timestep of 5 min, and, based on hypotheses we test later, only explicitly modelled the turning angles of the steps or ’bearing changes’. $\Delta \theta \left(t\right)$, the ‘bearing change’ at time t depends on the current bearing of the growth cone $\theta \left(t\right)$, the angle $\varphi \left(t\right)$ of the vector connecting the growth cone to its anchor point, the gradient direction Ψ and the noise ξ according to Equation (1):
where two parameters $a$ and $b$ scale the contributions of the first term representing persistence and the second term representing the bias due to the gradient. The symbol $\angle \left(x,y\right)$ denotes the angle difference $xy$ constrained to take values from −π to π. It is positive for an anticlockwise turn to get from $y$ to $x$. As the bearing is biased by the gradient direction, the overall growth cone angle $\varphi \left(t\right)$ will also be biased by the gradient, coupled through the above equation.
We first assume there is only one fixed anchor point where the axon initially grew out of the cell body or the main branch. We will later relax this assumption and allow the growth cone to put down new anchor points along its path. We denote the rate of anchor point deposition as r, which is the inverse of the average number of steps per new anchor point.
We first assume an initial direction of $\varphi \left(0\right)=\theta \left(0\right)=\pi \u22152$, a gradient direction of Ψ = 0, and a fixed step size s every 5 min. In the idealized noiseless case (ξ = 0) as $t\to \infty $, the equation reaches a steady state when $\Delta \theta =0$, that is:
This gives:
with $\alpha =\frac{a}{a+b}$. Defining L to be the distance of the growth cone from its original position, and using the geometry in Figure 1A, we have:
The approximation above is due to $s$ ≪ $L$ and ϕ(t) → Ψ = 0 as t → ∞. Using the Taylor expansion f(x_{0} + δx) ≈ f(x_{0}) + δxf’(x_{0}) and d tan^{−1}(x)/dx = 1/(x^{2} + 1), we invert both sides of the above equation to obtain:
At $t\to \infty $, $\Delta \theta \left(t\right)\to 0$, meaning the growth direction aligns with the gradient, thus $\varphi \left(t\right)\to 0$ and $L\approx st$ due to geometry (even for the $a$ = 0 case), so the above equation can be simplified as
Therefore, the longterm turning behaviour of axons in the model is given by the power law $\varphi \left(t\right)\propto {t}^{\left(\alpha 1\right)}$.
References

Quantifying mechanical force in axonal growth and guidanceFrontiers in Cellular Neuroscience 9:359.https://doi.org/10.3389/fncel.2015.00359

Growth cones as soft and weak force generatorsProceedings of the National Academy of Sciences of the United States of America 108:13420–13425.https://doi.org/10.1073/pnas.1106145108

The chemotactic effect of mixtures of antibody and antigen on polymorphonuclear leucocytesJournal of Experimental Medicine 115:453–466.https://doi.org/10.1084/jem.115.3.453

Branching patterns of individual sympathetic neurons in cultureThe Journal of Cell Biology 56:702–712.

Mechanical tension produced by nerve cells in tissue cultureJournal of Cell Science 37:391–410.

Axonal growth in response to experimentally applied mechanical tensionDevelopmental Biology 102:379–389.

Transport and turnover of microtubules in frog neurons depend on the pattern of axonal growthThe Journal of Neuroscience 18:821–829.

Molecular mechanisms of axon guidanceDevelopmental Biology 292:13–24.https://doi.org/10.1016/j.ydbio.2005.12.048

Microfluidic engineered high cell density threedimensional neural culturesJournal of Neural Engineering 4:159–172.https://doi.org/10.1088/17412560/4/2/015

Generation of gradients having complex shapes using microfluidic networksAnalytical Chemistry 73:1240–1246.

Molecular mechanisms of axon guidanceScience 298:1959–1964.https://doi.org/10.1126/science.1072165

GDNF acts as a chemoattractant to support ephrinAinduced repulsion of limb motor axonsCurrent Biology : CB 20:2150–2156.https://doi.org/10.1016/j.cub.2010.11.021

Genetic evidence for a contribution of EphA:ephrinA reverse signaling to motor axon guidanceThe Journal of Neuroscience 32:5209–5215.https://doi.org/10.1523/JNEUROSCI.570711.2012

Investigating axonal guidance with microdevicebased approachesThe Journal of Neuroscience 33:17647–17655.https://doi.org/10.1523/JNEUROSCI.327713.2013

The biophysics of neuronal growthReports on Progress in Physics 73:094601.https://doi.org/10.1088/00344885/73/9/094601

Autism spectrum disorders: developmental disconnection syndromesCurrent Opinion in Neurobiology 17:103–111.https://doi.org/10.1016/j.conb.2007.01.009

Cytomechanics of axonal developmentCell Biochemistry and Biophysics 27:135–155.https://doi.org/10.1007/BF02738107

Models of axon guidance and bundling during developmentProceedings of the Royal Society B: Biological Sciences 266:2231–2238.https://doi.org/10.1098/rspb.1999.0913

Signaling mechanisms in cortical axon growth, guidance, and branchingFrontiers in Neuroanatomy 5:62.https://doi.org/10.3389/fnana.2011.00062

Regulation of substrate adhesion dynamics during cell motilityThe International Journal of Biochemistry & Cell Biology 34:746–761.https://doi.org/10.1016/S13572725(01)001716

1433 proteins regulate protein kinase a activity to modulate growth cone turning responsesJournal of Neuroscience 30:14059–14067.https://doi.org/10.1523/JNEUROSCI.388310.2010

Mechanisms and molecules of neuronal wiring: a primerCold Spring Harbor Perspectives in Biology 3:a001727–14.https://doi.org/10.1101/cshperspect.a001727

A mathematical framework for modeling axon guidanceBulletin of Mathematical Biology 69:3–31.https://doi.org/10.1007/s1153800691424

Chemotactic response of nerve fiber elongation to nerve growth factorDevelopmental Biology 66:183–196.

BookAxonal pathfinding: Extracellular matrix roleIn: Spire LR, editors. Encyclopedia of Neuroscience, vol 1. Oxford: Academic Press. pp. 1139–1145.

Computer model of growth cone behavior and neuronal morphogenesisJournal of Theoretical Biology 174:381–389.

Asymmetric modulation of cytosolic cAMP activity induces growth cone turningThe Journal of Neuroscience 12:1253–1261.

Thalamocortical development: how are we going to get there?Nature Reviews. Neuroscience 4:276–289.https://doi.org/10.1038/nrn1075

The trip of the tip: understanding the growth cone machineryNature Reviews. Molecular Cell Biology 10:332–343.https://doi.org/10.1038/nrm2679

Growth factors and combinatorial therapies for CNS regenerationExperimental Neurology 209:313–320.https://doi.org/10.1016/j.expneurol.2007.08.004

New perspectives on neuronal development via microfluidic environmentsTrends in Neurosciences 35:752–761.https://doi.org/10.1016/j.tins.2012.09.001

Axon guidance by growthrate modulationProceedings of the National Academy of Sciences of the United States of America 107:5202–5207.https://doi.org/10.1073/pnas.0909254107

Time series characterization of simulated microtubule dynamics in the nerve growth coneAnnals of Biomedical Engineering 23:268–286.

Effect of nerve growth factor (nGF) on survival of superior cervical ganglion (sCG) transplanted into the third ventricle in rats japaneseJournal of Pharmacology 53:11–18.

Mathematical modeling of axonal formation. part i: geometryBulletin of Mathematical Biology 73:2837–2864.https://doi.org/10.1007/s1153801196482

Analysis of the growth cone turning assay for studying axon guidanceJournal of Neuroscience Methods 170:220–228.https://doi.org/10.1016/j.jneumeth.2008.01.014

Assays for eukaryotic cell chemotaxisCombinatorial Chemistry & High Throughput Screening 12:580–588.https://doi.org/10.2174/138620709788681952

Drosophila neurons actively regulate axonal tension in vivoBiophysical Journal 99:3208–3215.https://doi.org/10.1016/j.bpj.2010.09.029

Focal adhesion kinase signaling at sites of integrinmediated adhesion controls axon pathfindingNature Neuroscience 9:1274–1283.https://doi.org/10.1038/nn1762

A new chemotaxis assay shows the extreme sensitivity of axons to molecular gradientsNature Neuroscience 7:678–682.https://doi.org/10.1038/nn1259

A quantitative analysis of branching, growth cone turning, and directed growth in zebrafish retinotectal axon guidanceThe Journal of Comparative Neurology 521:1409–1429.https://doi.org/10.1002/cne.23248

What does the developing brain tell us about neural diseases?The European Journal of Neuroscience 35:1811–1817.https://doi.org/10.1111/j.14609568.2012.08171.x

The emerging role of forces in axonal elongationProgress in Neurobiology 94:91–101.https://doi.org/10.1016/j.pneurobio.2011.04.002

Cyclic nucleotidedependent switching of mammalian axon guidance depends on gradient steepnessMolecular and Cellular Neurosciences 47:45–52.https://doi.org/10.1016/j.mcn.2011.02.012

Colocalization of NGF binding sites, trk mRNA, and lowaffinity NGF receptor mRNA in primary sensory neurons: responses to injury and infusion of NGFThe Journal of Neuroscience 12:4011–4022.

A microfluidicsbased turning assay reveals complex growth cone responses to integrated gradientsof substratebound ECM molecules and diffusible guidance cuesLab Chip 586:227–237.

Neuronal and nonneuronal expression of neurotrophins and their receptors in sensory and sympathetic ganglia suggest new intercellular trophic interactionsThe Journal of Comparative Neurology 353:143–159.https://doi.org/10.1002/cne.903530113

Nerve growth cone guidance mediated by g protein–coupled receptorsNature Neuroscience 5:843–848.https://doi.org/10.1038/nn899

Navigating their way to the clinic: emerging roles for axon guidance molecules in neurological disorders and injuryDevelopmental Neurobiology 67:1216–1231.https://doi.org/10.1002/dneu.20512

Tensile regulation of axonal elongation and initiationThe Journal of Neuroscience 11:1117–1125.

Ability of polymorphonuclear leukocytes to orient in gradients of chemotactic factorsThe Journal of Cell Biology 75:606–616.
Article and author information
Author details
Funding
National Health and Medical Research Council (1083707)
 Geoffrey J Goodhill
Gatsby Charitable Foundation
 Peter Dayan
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Ethics
Animal experimentation: Experiments involving animals were conducted in accordance with the Animal Care and Protection Act Qld (2002) and the Australian Code of Practivce for the Care and Use of Animals for Scientific Purposes, 7th edition (2004). Ethics approval was obtained from The University of Queensland Ethics Committee, approval number AE02608.
Version history
 Received: October 12, 2015
 Accepted: December 18, 2015
 Version of Record published: February 2, 2016 (version 1)
 Version of Record updated: March 1, 2016 (version 2)
Copyright
© 2016, Nguyen et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 2,569
 views

 450
 downloads

 988
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Neuroscience
Pavlovian fear conditioning has been extensively used to study the behavioral and neural basis of defensive systems. In a typical procedure, a cue is paired with foot shock, and subsequent cue presentation elicits freezing, a behavior theoretically linked to predator detection. Studies have since shown a fear conditioned cue can elicit locomotion, a behavior that  in addition to jumping, and rearing  is theoretically linked to imminent or occurring predation. A criticism of studies observing fear conditioned cueelicited locomotion is that responding is nonassociative. We gave rats Pavlovian fear discrimination over a baseline of reward seeking. TTLtriggered cameras captured 5 behavior frames/s around cue presentation. Experiment 1 examined the emergence of dangerspecific behaviors over fear acquisition. Experiment 2 examined the expression of dangerspecific behaviors in fear extinction. In total, we scored 112,000 frames for nine discrete behavior categories. Temporal ethograms show that during acquisition, a fear conditioned cue suppresses reward seeking and elicits freezing, but also elicits locomotion, jumping, and rearing  all of which are maximal when foot shock is imminent. During extinction, a fear conditioned cue most prominently suppresses reward seeking, and elicits locomotion that is timed to shock delivery. The independent expression of these behaviors in both experiments reveal a fear conditioned cue to orchestrate a temporally organized suite of behaviors.