1 Introduction↩︎

The human brain comprises an immense network of neurons whose axons and dendrites (collectively known as neurites) establish long range, highly specific connections during development [1]–[5]. Axons often extend over distances of tens to hundreds of cell diameters to reach appropriate targets, a process directed by the growth cone, a dynamic sensor and actuator complex at the axon tip that integrates biochemical, mechanical, and topographical cues to drive directed motion [4], [6]–[9]. Guidance signals include diffusible molecules (e.g., Netrins, Slits, Semaphorins) and substrate bound factors (Ephrins, extracellular matrix components, adhesion molecules), together with physical inputs such as stiffness, geometry, and electric fields [3], [9]–[16]. Notably, growth cones navigate heterogeneous microenvironments with high precision, continuously probing and updating their trajectory [4], [15], [17], [18]. These processes ultimately wire circuits that enable reflexes, learning, attention, and memory.

A central challenge is to translate this complexity into predictive, quantitative laws. Axonal extension is inherently noisy: ligand detection near the single-molecule limit, stochastic reaction networks, intermittent adhesion engagement, and fluctuating cytoskeletal remodeling all introduce variability at the scale of the growth cone. At the same time, cells exploit regulatory feedback to stabilize motion and amplify relevant cues [18]–[22]. In the molecular clutch model, actin polymerization at the leading edge, retrograde flow driven by myosin-II, and dynamic coupling to substrate adhesions (integrins, cadherins) together determine traction forces and growth cone advance [3], [4], [6], [23]–[31]. Positive feedback reinforces forward motion and alignment, whereas negative feedback damps fluctuations and prevents uncontrolled responses [4], [11], [12], [18]. These considerations strongly motivate modeling frameworks based on stochastic processes, such as Langevin and Fokker–Planck equations, that explicitly couple deterministic axonal guidance (drift) with random fluctuations (diffusion), and inherently incorporate feedback-modulated noise. In these models, clutch‑mediated force transmission contributes to the drift term, while polymerization of the cytoskeleton, adhesion turnover, and signaling noise set the diffusion term.

Stochastic differential equations (SDE) provide a compact way to encode this interplay: drift fields represent guidance and feedback-regulated tendencies (e.g., alignment torques induced by patterned substrates), while diffusion coefficients capture intrinsic and extrinsic noise whose magnitude may depend on the local microenvironment. From a modeling perspective, axonal trajectories reflect the interplay of (i) deterministic biases imparted by chemical gradients, substrate mechanics, and geometry, and (ii) stochastic fluctuations arising from receptor binding, signaling cascades, adhesion engagement, and cytoskeletal remodeling. This formalism supports parameter inference from data and enables hypothesis tests about mechanism (e.g., clutch-mediated traction vs. gradient sensing) using likelihood-based or information-theoretic criteria [19]–[21], [32]–[37].

Recent advances in microfabrication and microfluidics provide controlled in vitro platforms to calibrate and test theoretical models. Engineered culture systems allow independent tuning of biochemical, mechanical, and geometric cues and reveal strong stiffness dependence of axonal elongation and robust guidance by patterned topographies [22], [27], [29], [30], [38]–[44]. For example, on grooved or anisotropic substrates, axons display biased, persistent motion and enhanced alignment (Figure 1 shows an example of alignment for axonal growth on micropatterned surfaces). These are signatures naturally captured by drift–diffusion, velocity-jump, or biased, persistent random walk models [22], [36], [45]–[51]. Such descriptions yield directly measurable predictions for speed and turning angle distributions, mean-square displacements, velocity and angular correlation functions, and motility coefficients, thereby linking single-cell biophysics (cytoskeletal dynamics, adhesion kinetics) to ensemble-level trajectory statistics and, ultimately, circuit-level connectivity.

In our previous studies, we have demonstrated that cortical neurons grown on poly-D-lysine (PDL) coated substrates with periodic micropatterns exhibit strong alignment with the surface features [45]–[49], [52]. This behavior is consistent with an effective, substrate-induced deterministic torque and with feedback that modulates adhesion and cytoskeletal dynamics [22], [31]. We quantified axonal speeds, angular distributions, autocorrelation functions, diffusion (cell motility) coefficients, and cell–surface interaction forces. We have also extracted mechanical parameters such as elastic and bending modulus relevant to shape control during guidance [22], [42], [49]. Datasets obtained from these experiments are well suited for building stochastic models, testing scaling predictions (e.g., mean squared displacement vs. time), and quantifying the dependence of drift and diffusion on controllable environmental cues.

A quantitative, stochastic account of axonal guidance is also essential for engineering growth-permissive microenvironments and for therapeutic strategies in nerve repair and neurodegeneration. Insight into how feedback and noise shape axonal outgrowth can inform the design of neuroprosthetic scaffolds and bioinspired platforms that promote targeted regeneration and functional reconnection [38], [44], [53]–[59]. Because stochastic frameworks return experimentally measurable parameters (drift strengths, diffusion coefficients, correlation times), they provide actionable targets for materials design (e.g., stiffness, biochemical composition, pattern periodicity) and for pharmacological modulation of adhesion and cytoskeletal dynamics. For example, increasing pattern anisotropy strengthens orientation drift, whereas stabilizing adhesions reduces effective diffusion by lengthening correlation times [22], [49], [50].

This paper advances a stochastic perspective on neuronal growth. We present evidence that mechanical and biochemical guidance cues operate through feedback-regulated mechanisms to produce biased yet noisy trajectories, and we formalize this behavior with Langevin and Fokker–Planck models. We show how to use these models to interpret experimental data and to extract diffusion coefficients, drift fields, turning-angle and speed distributions, and parameters quantifying cell–substrate coupling. Together, these results demonstrate that stochastic models are not merely descriptive: they provide a compact, predictive framework linking intracellular dynamics to axonal guidance and network formation in complex, fluctuating environments. We begin with a general overview of mathematical models of neuronal growth. We then formulate Langevin and Fokker–Planck models with nonlinear drift terms, analyze growth angle and velocity distributions, and calibrate parameters using trajectory data on micropatterned substrates. Finally, we introduce a stochastic mechanochemical model that couples a Fokker–Planck description of growth cone position to explicit dynamics of actin polymerization, myosin-II–driven contraction, and point-contact adhesions. We show that regulated feedback loops between polymerization, contractility, and adhesion govern growth velocity and adhesion-dependent traction, collectively producing steady extension, damped oscillations, and limit cycles. A linear stability analysis of the coarsed-grained dynamics delineates the parameter regimes associated with each behavior and clarifies how these coupled feedbacks jointly tune axonal outgrowth.

2 Mathematical Modeling of Axonal Growth↩︎

Modeling frameworks for axonal growth range from phenomenological descriptions of trajectories to mechanistic models of growth cone biophysics. Early work treated growth-cone motion as a random walk, asking whether observed elongation–retraction dynamics could be explained without invoking detailed intracellular mechanisms. For example, Katz and colleagues showed that net advance and retraction are, to first approximation, consistent with an uncorrelated random walk [60], whereas Odde and collaborators reported short-time correlations between extension and subsequent retraction on minute timescales [61]. In parallel, Buettner and co-workers extracted probabilistic rules for filopodial extension–retraction from time-lapse imaging and formalized these rules into a stochastic model [62], [63]. At the level of chemical sensing, the Goodhill group developed statistical models of cue–receptor binding at the growth cone [35], deriving constraints set by gradient shape and noise on detectability and steering [64], and showing that spatial sensing outperforms purely temporal strategies across experimentally relevant concentration ranges [65]. Katz and Lasek further identified constraints required to obtain ordered axonal ensembles from simple random-walk processes [66]. These studies established the utility of stochastic kinematic descriptions and sensing-theory bounds for interpreting growth trajectories.

Because biochemistry and mechanics are multiscale, mechanistic modeling has concentrated on tractable subsystems or controlled environments. Segev and Ben-Jacob modeled self-wiring in diffusing guidance fields and used graph-theoretic metrics to compare emergent networks with experiments [67]. Van Ooyen’s group simulated multiple axons navigating domains with overlapping guidance cues [68]. At the subcellular level, Mogilner and Rubenstein developed a mechanical theory of filopodial architecture to infer optimal length and stability [69]. Padmanabhan and Goodhill incorporated a molecular feedback loop in cytoskeletal control pathways that yields unimodal or bistable outgrowth depending on point contact adhesion assembly rates. Combined with a stochastic angular process, this produces a random walk with rest model in which advance and pausing reflect the state of the internal switch [70]. Reduced compartmental models have also been proposed to predict growth cone responses to externally imposed gradients [71]. Collectively, these contributions link intracellular regulation (adhesion, cytoskeletal turnover, signaling feedback) to mesoscopic motion under explicit biophysical assumptions.

A complementary line of work casts axonal guidance as stochastic transport governed by SDEs, with deterministic drift encoding biases from chemical gradients, substrate mechanics, or geometry, and diffusion capturing intrinsic and extrinsic fluctuations (receptor noise, reaction networks, adhesion engagement, and cytoskeletal remodeling). Simulating these SDEs yields probability densities over position and orientation, enabling direct comparison with ensemble statistics (speed and turning-angle distributions, mean-square displacement, correlation functions) and testable predictions for competing biophysical mechanisms. Using such approaches, Hentschel and van Ooyen reproduced axonal bundling, guidance, and subsequent debundling in combined attractant–repellent fields [72]. Maskery and Shinbrot used simulations based on SDE to estimate minimum detectable gradients under realistic noise [73]. Pearson and colleagues obtaining baseline trajectory geometries for growth in cue-free environments, [32]. Goodhill and collaborators coupled ligand binding with filopodial dynamics to generate guided trajectories in imposed gradients [74]. At the subcellular scale, Betz and co-workers used SDE analysis to lamellipodial fluctuations, showing that observed bimodality emerges from actin-driven bistability [75]. These studies illustrate how drift–diffusion models provide compact, data-driven links between microenvironmental statistics and growth cone kinematics.

Beyond single-axon descriptions, agent-based and network-level models incorporate interaction rules (for example, fasciculation/defasciculation, competition for cues) and domain topology. With local sensing and adhesion rules, simulations recover collective alignment, bundle formation, and target selection in patterned or heterogeneous landscapes [67], [68]. In such settings, stochasticity is not merely noise but a resource: fluctuations enable escape from local traps, exploration of alternative routes, and sensitivity to weak biases, while feedback modulates noise levels to stabilize chosen paths [19], [21], [76].

3 Langevin and Fokker-Planck Formalisms for Modeling Axonal Dynamics↩︎

Axonal growth arises from the interplay between deterministic and stochastic components of growth cone motility. Deterministic biases emerge, for example, from preferred orientations imposed by substrate geometry, whereas the stochastic contributions originate from cytoskeletal polymerization (actin and microtubules), intracellular signaling, detection of low-concentration cues, biochemical reactions, and the formation and turnover of lamellipodia and filopodia [1]–[8]. Because of this interplay, single-neuron trajectories are not deterministically predictable. However, ensemble behavior can be captured by probability densities governed by the associated SDEs. In particular, Langevin dynamics and the associated Fokker-Planck equation (FPE) provide a compact framework for modeling axonal dynamics as drift–diffusion processes that integrate guidance cues with noise [32], [72]–[75].

These stochastic models are particularly powerful when calibrated and validated against controlled in vitro measurements. By fitting SDE/FPE parameters to axonal trajectories on engineered substrates, one can extract effective drift fields (mechanical guidance strengths, alignment torques), diffusion (cell motility) coefficients, and correlation times, then test scaling laws such as mean-squared displacement (MSD) growth or velocity and angular correlations. For example, in our prior work, we have shown that cortical neurons grown on PDL-coated glass exhibited dynamics consistent with an effective V-shaped potential that regulates growth rates [45]. On ratchet-like, tilted-nanorod (nano-ppx) surfaces, axons aligned along a preferred direction due to a substrate-induced deterministic torque. We have measured angular distributions and drift–diffusion coefficients that quantify this bias [77], [78]. In a separate set of experiments, we showed that periodic geometrical patterns impart strong directional bias to axonal growth (Figure 1) [22], [31], [46]–[49], [52]. We have measured growth cone speeds, velocity autocorrelations, axonal orientation distributions, diffusion coefficients, and neuron-substrate traction forces. These examples show how stochastic transport models serve as a unifying language to compare disparate conditions (chemical gradients, mechanical stifness, geometrical anisotropy) and to map microenvironmental control parameters to observable path statistics.

3.1 Langevin equation for axonal growth. In previous work [50] we have shown that axonal dynamics on uniform glass surfaces is described by an Ornstein-Uhlenbeck (Brownian) process, defined by a linear Langevin equation for the velocity \(\Vec{V}\): \[\frac{d\vec{V}}{dt} = -\,\gamma_g\,\vec{V} \;+\; \vec{\Gamma}(t) \label{eq:A1}\tag{1}\] with constant damping \(\gamma_g\) and Gaussian white noise \(\vec{\Gamma}(t)\).

From Equation 1 we calculate the mean square axonal length and the velocity autocorrelation function [50]. By comparing the theoretical predictions with the experimentally measured distributions for these parameters we can extract the two fundamental parameters that characterize the motion of axons on glass surfaces: the diffusion coefficient \(D\) and the characteristic time for the exponential decay of the velocity autocorrelation function \(\tau_g=1/\gamma_g\). For cortical neurons grown on PDL coated glass these parameters are : \(D = (16 \pm 2)\,\mu\mathrm m^{2}\,\mathrm h^{-1}\) and \(\gamma_g = (0.1 \pm 0.05)\,\mathrm h^{-1}\) [50].

Axonal dynamics on micropatterned surfaces is described by a non-linear Langevin equation: \[\frac{d\vec{V}}{dt} \;\equiv\; \vec{a}(\vec{V},t) \;=\; \vec{a}_{\;d}(\vec{V},t) \;+\; \vec{\Gamma}(t), \label{eq:A2}\tag{2}\] where \(\vec{a}_d\) is the deterministic component and \(\vec{\Gamma}(t)\) the stochastic term. The acceleration of axons is decomposed into a component parallel to the direction of motion \(\vec{a}_{d, \parallel}(\vec{V},t)\), and a component perpendicular to this direction \(\vec{a}_{d, \perp}(\vec{V},t)\) (inset in Figure 1). A separate analysis of the two motions leads to the following non-linear Langevin equations for the two components of the acceleration [50]: \[a_{\parallel}(V,t) \equiv \frac{d \vec{V}_{\parallel}}{dt} = a_0 |\mathrm{sin}(\theta)| \;-\; \gamma_1\,V \;-\; \gamma_2\,V^2 \;+\; \Gamma_{\parallel}(t) \label{eq:A3}\tag{3}\] \[a_{\perp}(V,t) \equiv \frac{d \vec{V}_{\perp}}{dt} = a_1 \,\cos\theta \;+\; \Gamma_{\perp}(t) \label{eq:A4}\tag{4}\] Here, \(\theta\) is growth angle, \(V\) is the growth cone speed, and \(a_0,a_1,\gamma_1,\gamma_2\) are velocity-independent parameters that characterize axonal dynamics on substrates with periodic geometries. \(\Gamma_{\parallel},\Gamma_{\perp}\) are independent Gaussian white noises for parallel and perpendicular growth. We have shown that all these parameters are experimentally measurable [50].

Equations 3 4 show that the axonal dynamics on surfaces with periodic geometries is described by non-linear Langevin equations, involving quadratic velocity terms and non-zero coefficients for the angular orientation of the growing axon. There are some very important consequences for axonal growth that follow from this type of dynamics. In particular, Equations 3 4 show angular alignment of axonal growth on micropatterned surfaces. The magnitude of the perpendicular acceleration \(a_{\perp}\) has a maximum value when the direction of axonal growth is perpendicular to the surface pattern (i.e., for \(\theta=\pi/2, \;3\pi/2\) in Figures  1, and  2), and it equals zero when the axon grows along the pattern (\(\theta=0, \;\pi\)). This shows that the perpendicular component of acceleration \(a_{\perp}\) tends to align the growth cone along the direction of the pattern. The net effect is that of a deterministic torque (quantified by the parameter \(a_{1}\)) which rotates the growth cone towards the surface geometrical pattern. Figure 2 shows examples of angular (Figure 2a) and speed (Figure 2b) distributions for axonal growth on a surface with periodic micropatterns.

Another prediction of the model described by Equations 3 4 is that the growth cone reaches a terminal speed along the direction of the pattern, which can be found from the condition that the average acceleration in Equation 3 equals zero. This gives the following analytic expression for the terminal speed of the growth cone: \[V_{\rm ter} \;=\; \sqrt{\frac{a_0}{\gamma_2}\cdot|{\sin\theta}|+\frac{\gamma_1 ^2}{4\gamma_2 ^2}} -\frac{\gamma_1}{2\gamma_2} \label{eq:A5}\tag{5}\]

Equation 5 has a number of features that can be tested experimentally. First, the growth cones reaches terminal speed only for growth angles \(\theta \neq 0\). In addition, the terminal speed depends only on the ratios of the growth parameters \(a_o/\gamma_2\), and \(\gamma_1/\gamma_2\) which ultimately depend on the substrate geometry. Another important consequence of the non-linear Langevin Equations 3 4 is that axonal growth displays a cross-over from Brownian motion at earlier to a supper-diffusion regime at later times. The supper-diffusive dynamics is characterized by non-Gaussian speed distributions and power law increase of the axonal mean square length with time [49], [52]. From a biological perspective, the observed transition between the diffusive to super-diffusive axonal motion suggests long-range spatial and temporal correlations in the underlying dynamics [49].

3.2 Fokker-Planck equations for axonal growth. In a series of papers [22], [46], [49], [50] we have shown that the axonal dynamics on surfaces with periodic geometries is completely described by the following system of Fokker-Planck equations.

3.0.0.1 (a) Fokker–Planck equation for spatial probability \(P(\vec{r,t)}\) (Smoluchowski form):

\[\partial_t P(\vec{r},t) = D\,\nabla^2 P(\vec{r},t) \;+\; \frac{1}{\gamma}\,\nabla\!\cdot\!\bigl(P(\vec{r},t)\,\nabla U(\vec{r})\bigr) \label{eq:A6}\tag{6}\] with diffusion (motility) coefficient \(D\), damping \(\gamma\), and effective potential \(U(\vec{r})\). The 1D stationary solution along micropatterns is [22]: \[P(\vec{r},t)=A\,p(x,t),\qquad p(x)=\frac{1}{Z}\exp\!\left[-\frac{\gamma}{D}\,U_{\rm eff}(x)\right], \quad U_{\rm eff}=U_{\rm sub}(x)+U_{\rm fb}(x)+U_{\rm int}(x) \label{eq:A7}\tag{7}\] where \(A,Z\) are normalization constants, \(U_{\rm sub}\) is the external potential imposed by the substrate geometry, \(U_{\rm fb}\) the feedback potential, and \(U_{\rm int}\) accounts for neuron–neuron interactions. The form of these three potentials has been studied in reference [22].

3.0.0.2 (b) Fokker–Planck equation for the speed distribution \(P(v,t)\):

\[\partial_t P(v,t) = \partial_v\!\Big[\gamma_v\,(v-\bar v)\,P(v,t)\Big] \;+\; \frac{\sigma^2}{2}\,\partial_v^2 P(v,t) \label{eq:A8}\tag{8}\] The solution of this equation for the stationary speed distribution is: \[P(v)=B\,\exp\!\left[-\,\frac{\gamma_v}{\sigma^2}\,\bigl(v-\bar v\bigr)^2\right], \qquad \int_{0}^{\infty}P(v)\,dv=1 \label{eq:A9}\tag{9}\] with damping \(\gamma_v=1/\tau\) (relaxation time \(\tau\)), mean speed \(\bar v\), Gaussian noise strength \(\sigma\), and normalization constant \(B\).

3.0.0.3 (c) Fokker–Planck equation for the angular probability \(P(\theta, t)\):

\[\partial_t P(\theta,t) = \partial_\theta\!\Big[-\,\gamma_\theta\cdot\cos\theta\cdot P(\theta,t)\Big] \;+\; D_\theta\,\partial_\theta^2 P(\theta,t) \label{eq:A10}\tag{10}\] The solution of this equation for the stationary angular distribution is: \[P(\theta) = C\,\exp\!\Big[\frac{\gamma_{\theta}}{D_{\theta}}|\sin\theta|\Big], \quad \int_{0}^{2\pi}P(\theta)\,d\theta=1 \label{eq:A11}\tag{11}\] In Equations 10 11 \(P(\theta,t)\) is the probability distribution for the growth angle \(\theta\), \(C\) is a normalization constant, \(D_{\theta}\) represents the effective angular diffusion coefficient, and \(\gamma_{\theta}\cos \theta(t)\) corresponds to a “deterministic torque” representing the tendency of the growth cone to align with the preferred growth direction imposed by the surface geometry. The absolute value \(\lvert \sin\theta \rvert\) in Equation 11 reflects the symmetry of axonal growth around the \(x\) axis: the angular distributions centered at \(\theta\) are symmetric with respect to the directions \(\theta\) and \(\pi-\theta\) (Figure 2). This is a consequence of the fact that there is no preferred direction along the substrate micropattern. We also note that the deterministic torque has a maximum value if the growth cone moves perpendicular to the surface patterns (\(\theta=0\) or \(\theta=\pi\)), in which case the cell–surface interaction tends to align the axon with the surface pattern. The torque is zero for an axon moving along the micropattern.

In previous work, we have used the model given by Equations 6 11 to extract key dynamical parameters of axonal motion. Typical values obtained for the growth parameters are: diffusion coefficient \(D=(22\pm4)\,\mu\mathrm{m}^2/\mathrm{hr}\), coefficient for the “deterministic” alignment torque \(\gamma_{\theta}=(0.13\pm0.04)\,\mathrm{hr}^{-1}\), and characteristic time for axonal alignment \(\tau=(5.1\pm0.8)\,\mathrm{hr}\). We have performed a detailed analysis of how these parameters depend on the type of substrate, growth time, and chemical modification of the neurons [22], [46], [47]. These results show that the dynamics of the ensemble of axons can be described phenomenologically if each growth cone is modeled as an automatic controller with a closed feedback loop [22]. Growth alignment is fully determined by the surface geometry, and the distance between micropatterns plays the role of a control parameter. In particular, we have performed experiments which demonstrate that the disruption of cytoskeletal dynamics through neuronal treatment with different chemical compounds alters the feedback loop of the cellular controller [22], [46], [47].

4 Mechanical Beam Model of Axons↩︎

The phenomenological models discussed in the previous sections form a basis for quantifying cell-cell and cell-surface interactions, and ultimately for describing how the formation of neuronal network emerges from collective biophysical processes of single cells. In particular, the Fokker-Planck dynamics can be justified by a simple mechanical model that takes into account the cell-substrate interactions [22]. The model considers the bending-induced strained sustained by the axon while growing on the semi-cylindrical pattern of radius \(R\): axonal adhesion to the surface leads to axonal bending, which in turn leads to increased mechanical strain energy in the axon cytoskeleton. The mechanical strain energy \(E\) depends on the axon bending modulus \(F\), and the local surface curvature \(\kappa(\theta,R)\) [22] : \[E \;=\; \frac{F}{2}\,\kappa^2(\theta,R) \label{eq:A12}\tag{12}\] In the case of axonal growth on the micropatterned surfaces, the curvature of an axon around the cylindrical pattern of radius \(R\) is given by: \[\kappa(\theta,R) \;=\; \frac{|\cos\theta|}{R} \label{eq:A13}\tag{13}\] For the stationary growth described by the Fokker-Planck model one can assume a Boltzmann-type distribution for the probability of axon growing in a given direction: \[P(\theta)\;=\; A_1\,\exp\!\left(-\frac{E}{E_0}\right) \;= A_1\,\exp\!\left(-\frac{F}{E_0\cdot R^2}\cdot\cos^2\theta \right) \; \label{eq:A14}\tag{14}\] where \(E_0\) is the characteristic energy scale for axonal bending, and \(A_1\) is an overall normalization constant. By comparing the solution of this simple mechanical beam model (Equation 14 ) with the stationary solutions of the Fokker–Planck model (Equations 9 and 11 ), and using the experimentally measured values for the radius \(R\) of curvature of the micropattern and the growth parameters \(D_{\theta}\) and \(\gamma_{\theta}\), one can extract the bending modulus of the axon. Typical values for the bending modulus are \(F \approx 23~\mathrm{J}\,\mu\mathrm{m}^{2}\) for untreated neurons, and \(F \approx 17~\mathrm{J}\,\mu\mathrm{m}^{2}\) for neurons in which the cytoskeletal dynamics was inhibited by chemical modification [22].

These results indicate that axonal stiffness and substrate curvature can act jointly to direct axonal growth. The framework can be extended to include explicit dependencies of the growth parameters on biomechanical and geometric guidance cues, such as substrate geometry and stiffness, as well as on externally applied forces. For example, a proposed model for the cooperative motion of close‑packed cell ensembles treats contractile forces and effective cellular polarization as internal variables that generate waves of collective migration [79]. Continuum mechanical models that couple cell–substrate interactions to cellular biomechanical properties have likewise been proposed [80]–[83]. The parameters in these models are accessible experimentally via combined Atomic Force Microscopy (AFM) and Traction Force Microscopy (TFM) measurements [22], [31], [42].

5 Linear Stability and Oscillations in a Stochastic Actin-Myosin-Clutch Model↩︎

The preceding sections established a drift–diffusion description for growth cone motion using Langevin dynamics and the associated FPE, with drift fields encoding guidance and feedback, and diffusion terms capturing intrinsic and extrinsic noise. We now build on this framework by coupling the growth cone kinematics to a minimal mechanochemical model for actin polymerization, myosin–II contractility, and point–contact adhesions. This model links the probabilistic description of axon position to explicit intracellular processes, allowing us to analyze when feedback produces steady outgrowth, damped oscillations, or sustained limit cycles.

Experimental studies have established that cell-substrate adhesion is mediated by cell adhesion molecules (CAMs) which act as “molecular clutches” that couple the actin filaments (F-actin) to the underlying substrate. Growth cones form complex point contact adhesions, which are hierarchical assemblies where integrin molecules bind to extracellular matrix proteins on the growth substrate [2]–[4]. Simultaneously, numerous other CAMs, including those involved in signal transduction and actin binding, are recruited to the intracellular side [4]–[6]. According to the clutch model, myosin contracts actin filaments and generates active contractile stress, while adhesion complexes produce frictional adhesion forces opposing the movement of actin filaments [6], [11], [12]. We assume that the growth velocity \(v_g\) depends on three coarse–grained internal variables: polymerized actin density \(\rho_a(x,t)\), bound myosin density \(M(x,t)\), and point–contact adhesion density \(A(x,t)\).

Assuming that the contractile stress is proportional to the density \(M\) of myosin bound to F-actin, we use the following force balance equation for the local velocity of the actin flow in the lab coordinate system [83]: \[\nabla \cdot \left[ \mu_A \nabla v_{\text{actin}} + k M \right] = \xi(A, \rho_a) v_{\text{actin}}, \label{balance}\tag{15}\] The first term in Equation  15 is the sum of the divergences of the passive shear and deformation stresses for F-actin (where \(\mu_A\) denotes the actin viscosity). The term \(kM\) represents the active contractile force (with \(k\) the force per myosin unit), and \(\xi\!\left(A,\rho_N\right)\, v_{\textrm{actin}}\) is the frictional adhesion force. The coefficient \(\xi(A, \rho_a)\) measures the strength of growth-cone–substrate adhesion and depends on the density \(A\) of point-contact adhesions and on the actin density \(\rho_a(x,t)\). The exact functional dependence has not been characterized in the literature. Following earlier work on keratocytes [82], we model adhesion with a power law dependence: \(\xi(A, \rho_a) = \xi_0 A^\alpha \rho_a^\beta\), where \(\xi_0\) is the adhesion coefficient and \(\alpha\) and \(\beta\) are feedback exponents.

Consistent with the clutch models, the actin polymerization rate scales linearly with actin density, \(v_p = \nu \rho_a\). The net growth cone velocity is then given by [82], [83]: \[v_g = \frac{\nu \rho_a}{1 + \frac{\xi_0 A^\alpha \rho_a^\beta}{k}} -\frac{kM}{\xi_0 A^{\alpha}\rho_a^{\beta}} ,\] and the Fokker-Planck Equation 8 yields a Gaussian-like steady-state velocity distribution: \[P_{\text{s}}(v) \propto \exp\left[-\frac{(v - \langle v_g \rangle)^2}{2\sigma_v^2}\right],\] where \(\sigma_v^2 \sim D_g\) (effective diffusion coefficient) captures the stochastic fluctuations.

To perform the stability analysis, we adopt a coarse grained (spatially averaged) description for \((\rho_a,M,A)\) and retain the dominant kinetic couplings. Actin, myosin, and point contact dynamics are then described by [82], [83]: \[\frac{d\rho_a}{dt}= -\gamma \rho_a - v_g \frac{\partial \rho_a}{\partial x} \label{eq:rho}\tag{16}\]

\[\frac{dM}{dt} = -k_{\text{off}} M + k_{\text{on}} (M_0 - M) \label{eq:myo}\;\tag{17}\]

\[\frac{dA}{dt} = \alpha_A - \beta_A A \label{eq:adh}\tag{18}\] where \(\gamma\) is the rate of actin disassembly in the growth cone; \(k_{\textrm{off}},\;k_{\textrm{on}}\) denote the myosin detachment and attachment rates, respectively; and \(M_0\) is the concentration of unbound myosin (which attaches to F-actin). We also assume that point contact adhesion complexes appear across the growth cone with constant rate \(\alpha_A\) and disassemble with rate \(\beta_A\).

For stability analysis we linearize Equations 16 –18 around the steady state \((\rho_a^*, M^*, A^*)\). Defining the perturbations: \[\delta\rho_a = \rho_a - \rho_a^*, \quad \delta M = M - M^*, \quad \delta A = A - A^*\] the system  16 –18 can be written as: \[\frac{d}{dt} \begin{bmatrix} \delta \rho_a \;\delta M \;\delta A \end{bmatrix} = \mathcal{J} \begin{bmatrix} \delta \rho_a \;\delta M \;\delta A \end{bmatrix},\] where \(\mathcal{J}\) is the Jacobian matrix. These terms form a feedback loop that quantifies how changes in each internal variable modulate the axonal drift velocity \(v_g\). The local dynamics are determined by the eigenvalues of \(\mathcal{J}\), which satisfy a cubic characteristic equation: \[\lambda^{3}+a_{1}\lambda^{2}+a_{2}\lambda+a_{3}=0, \label{eq:char}\tag{19}\] Routh–Hurwitz stability criterion imply linear stability when \(a_{1}>0\), \(a_{2}>0\), \(a_{3}>0\), and \(a_{1}a_{2}>a_{3}\) [34], [84]. A Hopf bifurcation (the onset of small, self–sustained oscillations) occurs when \(a_{1}>0\), \(a_{2}>0\), \(a_{3}>0\) and \(a_{1}a_{2}=a_{3}\). In the present model, oscillations arise when clutch reinforcement (large \(\xi(A, \rho_a )\)), load–sensitive myosin recruitment (intermediate \(k_{\textrm{on}}\)), or strong adhesion/actin feedback exponents (\(\alpha,\beta\)) amplify the effective gain in the \(\rho_a\to v_g \to M,A\to v_g\) loop, so that the trace of \(\mathcal{J}\) remains negative but \(a_{1}a_{2}\) approaches \(a_{3}\). Conversely, rapid depolymerization (\(\gamma\) large), fast adhesion turnover (\(\beta_A\) large), or weak mechanosensitive feedback (small \(\xi_0\)) push the system deep into the stable regime with purely real, negative eigenvalues (overdamped return to steady growth). A summary of the dynamical behavior is provided in Table 1.

Illustrative eigenvalue spectrum. As a concrete example, we consider a reduced mechanochemical model of growth cone dynamics. Let \(x\) denote the dimensionless, actin-driven protrusion (tip) velocity, normalized by a characteristic polymerization speed. This quantity captures contributions from actin polymerization, myosin contractility, and cytoskeletal resistance. Let \(y\) denote the normalized adhesion density, proportional to the number of substrate-bound point contact adhesions (e.g., integrin complexes), and scaled by a reference density. The pair (\(x\), \(y\)) constitutes the minimal state vector for modeling the coupled dynamics of axonal extension and adhesion regulation under feedback control. For the dynamical variables \(x\) and \(y\), the Jacobian is: \[J = \begin{pmatrix} \frac{\partial \dot{x}}{\partial x} & \frac{\partial \dot{x}}{\partial y} \\ \frac{\partial \dot{y}}{\partial x} & \frac{\partial \dot{y}}{\partial y} \end{pmatrix}_{(x^*,y^*)},\] and we define the effective bifurcation parameter \(\mu\) as the trace of the Jacobian matrix of the linearized system at the steady state \((x^*,y^*)\) [34], [84]: \[\mu \;=\; \mathrm{Tr}\,J \;=\; \left.\frac{\partial \dot{x}}{\partial x}\right|_{(x^*,y^*)}+ \left.\frac{\partial \dot{y}}{\partial y}\right|_{(x^*,y^*)}.\]

In the reduced mechanochemical model of actin-driven protrusion and clutch-mediated force transmission, the dimensionless bifurcation parameter \(\mu\) can be interpreted as a composite control parameter encoding the balance between actin polymerization, substrate adhesion, and myosin-driven contractility [34], [76]. Specifically, \(\mu\) increases with higher actin polymerization rate, which drives forward motion; stronger effective adhesion \(\xi(A,\rho)\), which promotes traction force buildup; and reduced myosin contractility \(k \nabla M\), which otherwise resists forward motion. Thus, \(\mu\) acts as an effective propulsion-to-friction ratio. The phase portraits for the reduced mechaniochemical model are shown in (Figure 3). Negative values of \(\mu\) correspond to a stable focus with damped protrusion dynamics (Figure 3 (a)); \(\mu=0\) marks a Hopf bifurcation (Figure 3 (b)); and positive \(\mu\) indicates an unstable focus or limit-cycle oscillations (Figure 3 (c)).

Figure 4 shows the corresponding temporal profiles of the dimensionless growth cone velocity \(v_g(t)\) for the same regimes of the control parameter \(\mu\) as in Figure 3. For \(\mu<0\), trajectories that spiral toward the fixed point in Figure 3 (a) manifest as monotonic relaxation or damped oscillations in time (Figures 4 (a) and (b)). Near the Hopf threshold (\(\mu\approx 0\)), the approach to a closed orbit in Figure 3 (b) becomes a stable limit cycle with sustained oscillations of \(v_g(t)\) (Figure 4 (c)). For \(\mu>0\), the outward spirals or repelling fixed points in Figure 3 (c) are reflected by exponential growth of \(v_g(t)\) (Figure 4 (d)). Together, these time traces confirm the bifurcation structure inferred from the eigenvalue analysis and provide observable readouts that distinguish stable, weakly underdamped, limit-cycle, and unstable regimes.

6 Discussion and Conclusions↩︎

The results presented here connect a stochastic description of axonal motion with an explicit mechanochemical account of growth cone regulation. In the first part of the paper, Langevin dynamics and the associated Fokker–Planck equation (FPE) provide a drift–diffusion framework for axonal trajectories, with drift fields encoding guidance and feedback and diffusion terms capturing intrinsic and extrinsic fluctuations. We then coupled this FPE drift to a minimal actin–myosin–clutch model that links polymerization-driven protrusion, myosin-II contractility, and point-contact adhesion. This integration yields a compact model in which measurable intracellular processes set the parameters of a stochastic transport equation, enabling direct comparison to trajectory statistics.

The linear stability analysis clarifies how positive and negative feedback loops jointly regulate robust yet adaptable outgrowth. Near a stable fixed point, linearizing the FPE with approximately constant effective diffusivity reduces the growth cone dynamics to an Ornstein–Uhlenbeck process. This predicts a Gaussian stationary distribution of growth cone velocities narrowly peaked around \(\langle v_g \rangle\) with variance controlled by effective diffusion coefficient \(D_g\). As parameters approach a Hopf bifurcation, the dominant eigenvalues of the coarse grained mechanochemical dynamics form a weakly damped complex pair. In this regime, stochastic trajectories exhibit quasi-oscillatory fluctuations, providing concrete experimental signatures of clutch-mediated feedback. Beyond the Hopf threshold, the model admits a stable limit cycle in the coarse-grained variables, consistent with self-sustained growth–retraction oscillations whose amplitude and frequency are set by the internal feedbacks rather than by initial conditions.

Biologically, these regimes map onto experimentally accessible control parameters. Strengthening adhesion reinforcement or increasing load-sensitive myosin recruitment raises the closed-loop gain, pushing the system toward oscillatory dynamics. In contrast, rapid actin turnover or fast adhesion disassembly lowers the gain and stabilizes steady extension. At the level of substrate mechanics and adhesion, increased anisotropy or substrate curvature that enhances alignment increases the drift magnitude and reduces angular diffusion, whereas softening the growth substrate or weakening integrin engagement increases effective diffusion and can suppress oscillations by limiting traction. This framework captures how coupled positive and negative feedback produce reliable morphogenic outcomes: positive feedback amplifies weak guidance cues into directed motion, while negative feedback bounds fluctuations and enables rapid adaptation to perturbations in biochemical and mechanical environments.

The model yields several testable predictions that link intracellular control to ensemble statistics. In the stable regime, the Ornstein-Uhlenbeck reduction implies exponentially decaying velocity autocorrelations with a single characteristic time and a Gaussian-like distribution of short-time displacements [52]. As the system approaches a Hopf bifurcation, the dynamics becomes oscillatory and the power spectrum acquires a narrow peak that sharpens as damping decreases. Pharmacological or genetic perturbations provide targeted validation: partial inhibition of myosin-II (reducing the contractile coefficient) should shift the spectrum toward lower frequency and diminish oscillation amplitude; interventions that enhance adhesion reinforcement (e.g., integrin activation) should have the opposite effect. The parameters required by the model—including effective motility coefficients, correlation times, and measures of cell–substrate coupling—can be estimated by combining trajectory analysis on patterned substrates with independent mechanical measurements (for example, atomic force and traction force microscopy measurements [31], [42]), enabling a quantitative model-experiment comparison.

Several limitations suggest directions for refinement. The coarse-grained mechanochemical variables capture dominant feedbacks but neglect spatial heterogeneity across the growth cone, time delays in signaling, and curvature-dependent kinematics. Extending the internal dynamics to spatially distributed fields would capture competition between local protrusion and retrograde actin flow and enable mode-selection analyses of spatiotemporal patterns. Likewise, the present FPE treats drift and diffusion as functions of coarse variables. When diffusion depends on position or orientation (multiplicative noise), the stochastic calculus convention (Itô versus Stratonovich) and any induced noise-interpretation drift (spurious-drift terms) should be specified and tested. At the population level, coupling single-axon dynamics through interaction rules (fasciculation/defasciculation, competition for cues) would bridge single-cell biophysics to the formation of axon bundles and network architecture. Finally, parameter identifiability warrants careful analysis: joint fits to displacement distributions, autocorrelations, and power spectra, together with independent mechanical measurements, provide complementary constraints that improve robustness and reduce model degeneracy.

In conclusion, by embedding a clutch-regulated mechanochemical model into the drift term of the Fokker-Planck equation, we obtain a predictive, general framework that links intracellular regulation to stochastic growth cone motion and observed trajectory statistics. The stability analysis identifies biologically meaningful regimes: steady extension, damped oscillations, and sustained limit cycles, and specifies how polymerization, contractility, and adhesion feedback tune transitions among them. Beyond clarifying mechanisms of axonal guidance and connectivity, this framework provides practical design rules for engineered substrates and neuroprosthetic scaffolds: pattern anisotropy and curvature modulate directional drift; adhesion chemistry and stiffness tune diffusion and closed-loop gain; and targeted perturbations shift the system toward or away from oscillatory regimes associated with exploratory growth. From an applied-mathematics standpoint, the framework unifies SDE/FPE modeling with linear and nonlinear stability analysis, and provides concrete routes for parameter estimation and control-oriented design. Together, these insights advance a quantitative framework for controlling axonal outgrowth in vitro and, ultimately, for promoting functional regeneration in engineered microenvironments.

acknowledgments↩︎

This research was funded by a Tufts Faculty Research Award (FRAC).

References↩︎