The multiscale mechanics of axon durotaxis

The mechanical interaction between the substrate and the cytoskeleton is mediated by adhesion proteins, henceforth generically termed cross-bridges. In the growth cone, tension is generated by the movement of actin which stretches these cross-bridges. This tension opposes the retrograde motion of actin and pulls on the substrate, resulting reciprocally in a reaction force acting upon the axon shaft. The cross-bridges randomly bind and unbind, with a rate that depends on their level of tension (Bell et al., 1984; Chan and Odde, 2008). Overall, random detachment and attachment of cross-bridges, combined with stretch-induced tension result in emergent frictional interactions between the actin and the substrate which depends on substrate stiffness (Chan and Odde, 2008; Sens, 2013). The emergence of this frictional force is the central mechanism considered here underlying growth cone locomotion.

To model the mechanical interaction between a single actin filament and the substrate, we adapt the approach of Sens (2013). The filament is modelled as a straight, infinite rigid rod moving relative to the substrate with velocity $V$, as shown in Fig. 3 (we also introduce the longitudinal position $X=Vt$). In the steady regime, $V$ is constant in time. We model the cross-bridges as linearly elastic springs connecting the actin filament to the substrate. Each spring has an individual tension $f$, elongation $x$, and spring constant $\kappa$ related through Hooke’s law by

The cross-bridge tension can be related to the time $\Delta t$ elapsed since its attachment. Indeed, for a perfectly rigid substrate, the tension of a single cross-bridge since attachment is $f\approx\kappa V\Delta t$ (assuming that cross-bridges form unstretched, and neglecting the geometric nonlinearity due to the distance between the actin filament and the substrate). The situation is more complex in the case of a deformable substrate, where the cross-bridges may interact indirectly and non-locally with one another through short- and long-range deformations of the substrate (Nicolas et al., 2008). This deformation may be expressed in terms of the Green solution for an infinite half-plane subject to surface traction (Hwu, 2010; Sens, 2013; Nicolas et al., 2008).

Here, for simplicity, we assume that the substrate is soft so that cross-bridges interact mechanically as the substrate deforms almost uniformly. The substrate is then modelled as a single spring $\kappa^{\prime}$ with which all cross-bridges are connected in parallel. To model the kinetics of cross-bridge detachment, we posit a Bell-type law (Bell, 1978), i.e. we assume that the probability that a cross-bridge detaches in a small time-interval $\Delta t$ is related to its tension $f$ according to

where $f_{0}$ is the characteristic force scale of detachment and $k_{\mathrm{off}}$ is the kinetic off-rate at rest. Further, we assume that cross-bridges form in the unstretched state at a constant rate $k_{\mathrm{on}}$, so the probability of attachment during $\Delta t$ is

Finally, we focus on a small section of the filament at the growth cone edge. Therefore, we assume that all material parameters are constant in space which reduces the problem to a zero-dimensional problem.

As can be seen, for a range of parameters, the traction-stiffness function $T\left({\chi}\right)$ may exhibit global maximum at a stiffness $\chi^{*}$. We call $\chi^{*}$ the optimal stiffness defined as the stiffness where traction is maximal $T_{\text{max}}=T\left({\chi}\right)$. At that value, the ability of the axon to develop traction forces on the substrate is maximal. From (21) we see that $T\left({\chi}\right)$ is maximal when $\mathcal{H}\left({u\left({\chi}\right)}\right)$ is maximal. The function $\mathcal{H}\left({x}\right)$ has a global maximum at $x=x^{*}\approx 0.387$ (at which point $\mathcal{H}(x^{*})\approx 0.298$), and therefore $T\left({\chi}\right)$ is maximised when $u\left({\chi}\right)=x^{*}$.

Substituting $u=x^{*}$ and $N=N^{*}$ into (20) and rearranging the terms, we obtain

This defines a system of two equations for the two unknowns $\left({\chi^{*},N^{*}}\right)$. Remarkably, there exists a closed-form solution for this system given by

with $\Gamma\left({s,x}\right)$ the upper incomplete gamma function.

Since $N^{*}$ is positive if $\chi^{*}$ is positive, we see from (23) that the condition for the existence of a positive local maximum $\chi^{*}$ is

When this condition is not satisfied, the traction is a strictly increasing function of $\chi$ with an horizontal asymptote as $\chi\rightarrow\infty$. Importantly, we will show in Section 3 that the emergence of negative durotaxis is directly predicated upon the existence of a positive optimal stiffness via the condition (24). Moreover, the second derivative $T^{\prime\prime}(\chi^{*})$ of the traction $T(\chi)$ at the optimal stiffness $\chi^{*}$ characterises the behaviour of the system near this optimum. Differentiating (21) twice with respect to $\chi$ and setting $\chi=\chi^{*}$ we obtain

An exact expression for $u^{\prime}\left({\chi^{*}}\right)$ in the previous equation is obtained by differentiating (20) w.r.t. $\chi$ at $\chi=\chi^{*}$, providing a system of two equations for the unknowns $u^{\prime}\left({\chi^{*}}\right)$ and $N^{\prime}\left({\chi^{*}}\right)$, which can be solved to obtain

with $\operatorname{Chi}$ and $\operatorname{Shi}$ the hyperbolic cosine and sine integrals, respectively.

The parameter $M$ quantifies the force of myosin contraction, which plays an important role in axon guidance (McCormick and Gupton, 2020; Turney and Bridgman, 2005). The parameter $M$ changes qualitatively the traction-stiffness curve. Indeed, we define $M_{\text{crit}}$ the value of $M$ at which the inequality (24) is an equality, that is,

This critical force defines two distinct behaviours.

First, for $M>M_{\text{crit}}$, when myosin contraction is sufficiently large, as shown in Fig. 4, there exists an optimal stiffness and the behaviour of the system is controlled by three important features of the traction-stiffness curve: (i) the optimal stiffness $\chi^{*}$; (ii) the curvature $T^{\prime\prime}\left({\chi^{*}}\right)$ defining how quickly traction changes close to the optimal stiffness; (iii) the asymptotic traction generated on a rigid substrate $T_{\infty}\vcentcolon=\lim_{\chi\rightarrow\infty}T\left({\chi}\right)$. We note that the optimal stiffness $\chi^{*}$ is a decreasing function of $M$ and tends to zero as $M\rightarrow\infty$; see (23). As can be seen from (25, 26), $T^{\prime\prime}\left({\chi^{*}}\right)$ is a bounded, increasing function of $M$ such that

This asymptotic value is the maximal curvature achievable at $\chi^{*}$ by modulating contractility $M$. We see that the curvature increases as $\chi^{*}\to 0$ as $M\to\infty$. However, there is a physical bound to the value of $M$ depending on available binding sites. This relationship shows that the larger the myosin contraction, the faster an axon will react to changes in stiffness close to the optimal stiffness.

To investigate the behaviour of $T_{\infty}$ in the limit $M\rightarrow\infty$, we first divide (20a) by $\chi$ and take the limit $\chi\rightarrow\infty$ to obtain

Note that $T_{\infty}=f_{0}K\mathcal{H}\left({u}\right)$ in the above equation and therefore we are concerned with the asymptotics of $\mathcal{H}\left({u}\right)$ in the limit $M\rightarrow\infty$. Showing that $u\sim\beta/M$ as $M\rightarrow\infty$, it follows that $T_{\infty}\sim f_{0}K\mathcal{H}\left({{\beta}/{M}}\right)\rightarrow 0$.

Second, when myosin contraction is weak, that is for $M<M_{\text{crit}}$, there is no maximum and traction is a monotonically increasing function of $\chi$, plateauing to $T_{\infty}=\lim_{\chi\rightarrow\infty}T\left({\chi}\right)$. Noting that $u\sim\left({\beta+f_{0}K}\right)/M$ as $M\rightarrow 0$ we obtain

We see that as $M$ increase past $M_{\text{crit}}$, $T_{\infty}$ decreases. Hence, for large myosin contraction, axons deposited on a very rigid substrate will exert a reduced traction.

Initially, a morphoelastic rod is defined by a smooth curve $\boldsymbol{R}_{0}\left({S_{0}}\right)=\left({X_{0}\left({S_{0}}\right),Y_{0}\left({S_{0}}\right)}\right)\in\mathbb{R}^{2}$ with arclength $S_{0}\in\left[0,L_{0}\right]$. We choose $S_{0}=0$ to be at the base (the soma) and $S_{0}=L_{0}$ at the axon tip (growth cone). As the axon grows and deforms elastically under loads, the rod adopts a new current configuration parameterised by the material coordinate $S_{0}$ and time $t$ by $\boldsymbol{r}\left({S_{0},t}\right)=\left({x\left({S_{0},t}\right),y\left({S_{0},t}\right)}\right)$. We define $s\left({S_{0},t}\right)$ the arclength in the current configuration, i.e. the contour length between the base and the material point that was at position $S_{0}$ and time $t=0$. The Frenet-Serret frame of the current curve is composed of the unit tangent vector $\boldsymbol{\tau}={\partial{\boldsymbol{r}}}/{\partial{s}}$ (pointing towards the growth cone) and the unit normal vector $\boldsymbol{\nu}$, perpendicular to $\boldsymbol{\tau}$. Both can be encoded by a single angle $\theta$ between the $x$-axis and the tangent vector so that

At the tip, the orientation of the axon agrees with that of the growth cone so that $\boldsymbol{\tau}(L_{0})=\boldsymbol{t}$ and $\boldsymbol{\nu}(L_{0})=\boldsymbol{n}$. Then the Frenet-Serret equations can be written as

where $\lambda$ is the stretch and $\kappa$ the signed curvature. The pair of scalar fields $\left({\lambda,\kappa}\right)$ then fully describes the geometry of the rod (stretching and bending).

Neglecting inertial effects and body torques, the balance of linear and angular momenta is expressed in terms of the internal resultant force $\boldsymbol{n}=(n_{x},n_{y})$ and moment $m$ applied by a cross-section $S_{0}^{+}$ onto a cross-section $S_{0}^{-}$ as (see chapter 5 of Goriely, 2017, for details)

where the apostrophe denoting differentiation w.r.t. $S_{0}$; and $\boldsymbol{b}=(b_{x},b_{y})$ is a density of body force per unit current length, that is used to model the adhesion between the growing shaft and the substrate as a linear frictional force with coefficient $\zeta$ (O’Toole et al., 2008; Oliveri et al., 2021; Oliveri and Goriely, 2022)

Here the overdot denotes the material derivative, linked to the Eulerian derivative through the formula $\boldsymbol{\dot{r}}-\dot{s}\boldsymbol{\tau}={\partial{\boldsymbol{r}(s,t)}}/{\partial{t}}$.

To close the system, we must formulate constitutive relationships that relate stresses to deformations. We assume that the rod is extensible but unshearable. In a growing rod, only part of the deformation results in stresses, reflecting the remodelling and growth of the material and the resulting change in reference configuration. To model this effect, we introduce an intermediate reference configuration $\boldsymbol{R}\left({S_{0},t}\right)$ with arclength $S(S_{0},t)$, representing the equilibrium configuration of the rod achieved upon removal of the loads. We decompose the total deformation multiplicatively, in terms of the elastic, and growth stretches, respectively $\alpha\vcentcolon={\partial{s}}/{\partial{S}}$ and $\gamma\vcentcolon={\partial{S}}/{\partial{S_{0}}}$:

Constitutively, we assume that only $\alpha$ and $\kappa$ generate stresses. Therefore, we posit a quadratic elastic energy density (energy per unit intermediate length) for a morphoelastic rod with intrinsic curvature $\hat{\kappa}$, of the form

where $K$ and $B$ are the extensional and bending stiffnesses, respectively. The corresponding constitutive relations are

with $n\vcentcolon=\boldsymbol{n}\cdot\boldsymbol{\tau}$, the longitudinal tension.

There are two remaining important components to model, the growth of the axonal shaft specified by $\gamma$, and the remodelling of the curvature $\hat{\kappa}$. Physically, growth is due to mechanical stresses along the shaft that cause the axon to creep irreversibly above some critical tension (Franze and Guck, 2010; Suter and Miller, 2011; Goriely et al., 2015), while new cytoplasmic material is incorporated simultaneously into the shaft (O’Toole and Miller, 2011; Oliveri et al., 2022). During this process, cross-sectional area remains mostly constant, indicating that the cell maintains a homeostatic density (Dennerll et al., 1989; Lamoureux et al., 1989); thus we assume that $K$ and $B$ are constant in time and uniform along the rod. To model diffuse growth at each point, we use the exponential growth law (Oliveri and Goriely, 2022)

where $\left({n-n_{0}}\right)_{+}=n-n_{0}$ for $n\geq n_{0}$ ($0$ otherwise); $n_{0}\geq 0$ is the critical tension below which no growth occurs ; and $\eta$ plays the role of a morphoelastic viscosity (in unit force times time). Henceforth, for simplicity, we neglect the threshold effect and we set $n_{0}=0$.

We also assume a relaxation process for the bending moment due to internal remodelling by using the curvature evolution law

with a characteristic timescale ${\tau_{\gamma}}$ of the same order as $\eta/n_{0}$.

The connection to the mesoscopic model is obtained through the boundary conditions. At the base $S_{0}=0$, the rod is clamped

where $\boldsymbol{r}_{0}$ and $\theta_{0}$ are the base position and angle. At its tip, the rod is subject to the traction exerted by the growth cone:

where $\boldsymbol{F}$ depends on the growth cone position and orientation through the mesoscopic law (32).

There are three important parameters in the problem, the dimensionless parameter

which measures the relative effect of adhesion w.r.t. viscosity, and the two characteristic lengths

The growth length $\Lambda$ measures the distance over which axonal tension and energy dissipate into the substrate due to friction (O’Toole et al., 2008; Oliveri et al., 2021), and $\mathcal{L}$ is the bending length that measures the typical curvature radius of a beam due to the typical transverse traction $F$ from the growth cone.

At any time, given $\gamma(S_{0},t)$ and $\hat{\kappa}(S_{0},t)$, the boundary-value problem formed by (36, 37, 38, 39, 42, 41) with the boundary conditions (45, 44), can be solved numerically to obtain the axon’s trajectory. However, the frictional interaction between the axon and its substrate dissipates traction forces over a region close to the growth cone which becomes short, relatively to the entire axon, as time progresses. Thus, only a short and essentially constant distal section of the axon is dynamic, rapidly causing important computational difficulties due to the Lagrangian specification of our problem. Rather than solving for the system directly, we can take advantage of this effect and focus our attention on the growth process localised in a distal region, and freeze the most proximal part of the axon to account for permanent adhesion between the axon and the substrate (and discard that frozen region from the problem). This process is described in the next section.

The existence of a characteristic lengthscale $\Lambda$ can be used to further simplify the problem. For axons, we have $\Lambda\ll\mathcal{L}$ as the force applied by the growth cone only creates small deflections (hence a large radius of curvature $\mathcal{L}$). This scaling has two important consequences.

First, the effect of the force applied at the tip of the axon by the growth cone on the rest of the axon vanishes quickly away from the tip. It implies that only a short zone at the tip of length $L$ that is of order $\Lambda$ is subject to a force. In terms of modelling, we exploit this result computationally by constraining the dynamic region of the axon to a moving active portion of fixed arclength $L$ while segments of the axon found at a greater arclength are fixed and do not need to be updated. In morphoelasticity, this process is referred as tip growth (Goriely, 2017) and is found in plant shoots or roots, which undergo elongation over an effectively finite apical (or subapical) region (Silk and Erickson, 1979).

Second, since the active zone is small with respect to the typical radius of curvature, only small deflections away from the natural state take place in response to the growth cone. The distal end of the beam–the axon’s tip–is subject to the locomotory force $\boldsymbol{F}$ produced by the growth cone. The basal end makes the junction between the dynamic region and the frozen region, and is assumed to be clamped from the point of view of the elastic quasi-equilibrium. Therefore, we can approximate the shape of the active zone by

where $\boldsymbol{\hat{r}}=\left({\hat{x},\hat{y}}\right)$ and $\boldsymbol{\hat{\tau}}(t)$ is the position of the clamp; $\boldsymbol{\hat{\tau}}=\cos\hat{\theta}\,\boldsymbol{e}_{x}+\sin\hat{\theta}\,\boldsymbol{e}_{y}$ and $\boldsymbol{\hat{\nu}}=-\sin\hat{\theta}\,\boldsymbol{e}_{x}\,+\cos\hat{\theta}\,\boldsymbol{e}_{y}$ are the tangent and normal vectors to the curve at the base, with $\hat{\theta}$ the angle of the clamp to the horizontal $\boldsymbol{e}_{x}$ (henceforth, we denote quantities evaluated at the base of the beam with a circumflex); $\sf{x}$ and $\sf{y}$, typeset in the sans-serif font, denote the coordinates of the rod in the frame of the clamp. For small deflection we have ${\sf x}(\sigma,t)\approx\sigma$, and the rod equations of the previous section simplify to the beam equations (Goriely, 2017, p. 147).

There is however an important subtlety appearing in the problem. Since the active zone is moving in time, we have to properly advect the position of the beam element as the axon grows. We keep track of the position and direction of the clamp through the position vector $\boldsymbol{\hat{r}}(t)$ and model the evolution of $\boldsymbol{r}(\sigma,t)$ as the limit in small time step $\Delta t$ of the following process (Fig. 6): At $t$ the beam is clamped at $\boldsymbol{\hat{r}}(t)$ with direction $\boldsymbol{\hat{\tau}}(t)$ and is subject to a distal load $\boldsymbol{F}(t)$ and viscous friction along its length. First, ${\sf y}(\sigma,t)$ evolves according to the linear beam equations for a time interval of length $\Delta t\ll 1$. Then, a linear segment of length $V(t)\,\Delta t$ is added to the end of the beam in the direction of the tip. We call $V(t)$ the growth rate of the axon. Finally, the length constraint is enforced by setting the new profile of the beam, $\boldsymbol{r}(\sigma,t+\Delta t)$ to be the segment of length $L$ from the tip of the grown, deflected beam profile.

Taking the limit $\Delta t\rightarrow 0$, we obtain the continuous-time evolution equations (see details in A):

The deflection ${\sf y}(\sigma,t)$ satisfies the boundary conditions

Conditions (50a) define a clamp at $\sigma=0$, while (50b) enforces the load $\boldsymbol{F}$ applied with no torque on the beam at $\sigma=L$.

The growth speed $V$ depends on the force exerted by the growth cone on the tip of the axon, and is given, for small deflection and to leading order, by (O’Toole et al., 2008; Oliveri and Goriely, 2022; Oliveri et al., 2021)

We have developed a minimal model with two main components. From the microscopic model we have a mapping from substrate rigidity to the traction exerted by the growth cone, and from the mesoscopic model we obtain the shape of the axon as a function of time by locally integrating a clamped beam equation and advecting the clamp accordingly. We are now in a position to use this model to understand the effects of varying substrate rigidity on the motion and shape of axons.

We consider a heterogeneous field composed of two regions with different but uniform stiffnesses separated by a straight interface. In this case, the axon can either grow through the interface (refraction) or bounce against the interface (reflection), in analogy with optic ray theory (Oliveri et al., 2021). As the growth cone touches the interface, the force density along the growth cone margin becomes heterogeneous. The total force depends nonlinearly on the position and orientation of the growth cone via (31).

We assume that the axon grows from the region $x<0$, as shown in Fig. 7(a), and approaches the boundary with an angle of incidence $\theta_{1}\in[0,{\pi}/{2})$ (the angle between the axon and the normal to the interface). We assume that the interface is sharp but numerically, we approximate this straight interface by

where $T_{1}=T\left({\chi_{1}}\right)$ and $T_{2}=T\left({\chi_{2}}\right)$ are the force densities produced by the growth cone in the left-most and right-most regions, which depend on the stiffness of each region ($\chi_{1}$, $\chi_{2}$). The parameter $\lambda$ determines the length scale of transition from $T_{1}$ to $T_{2}$ across the interface and is chosen small enough with respect to the growth cone size that the results are not affected by further decreasing it (in practice, we found that $\lambda=R/4$ is suitable).

Initially, we assume a clamp position $\boldsymbol{r}_{0}\left({0}\right)$ far from the boundary in the incidence direction $\theta_{1}$. As soon as the growth cone leaves the interface region, it follows a straight trajectory and as $t\rightarrow\infty$, we define the refraction angle $\theta_{2}$ with the normal pointing into the right region. The axon may also reflect against the interface and stay in the left-most region as seen in Fig. 7(b). In this case, the reflection angle $\theta_{2}$ is defined as the angle between the left normal and the direction of the axon. Let $\theta_{\mathrm{m}}=\lim_{t\rightarrow\infty}\theta_{0}\left({t}\right)\mod 2\pi$, then restricting our attention to $\theta_{\mathrm{m}}\in[0,\pi)$ we have

By choosing the characteristic force-scale of the system to be $2T_{1}\sin{\varphi_{0}}$, which is equal to the magnitude of the total growth cone force when $\chi=\chi_{1}$, we obtain the two dimensionless parameters

In particular, $\beta_{a}$ relates to the typical radius of curvature of the axon, that is, smaller $\beta_{a}$ generates larger deflections.

We distinguish the two cases $\alpha_{T}<1$ and $\alpha_{T}>1$. When $\alpha_{T}>1$–see Fig. 7(a)–the durotactic force in the right region is greater than in the left region, thus the durotactic force component tends to pull the axon towards the right region. This results in a refractive behaviour with $\theta_{2}$ smaller than $\theta_{1}$.

At a superficial level, these observations are reminiscent of the refractive and reflective behaviours of light rays crossing a sharp-interface between two regions of different refractive indices. However, due to the finite size of the growth cone, the interaction with the interface is more complex and there is no simple quantitative analogy with Snell’s law. Further, a fundamental difference between this model and geometric optics is that we do not have time-reversibility in our system. In optics, the time-reversal symmetry of Maxwell’s equations implies that reflection occurs with $\theta_{1}=\theta_{2}$. As can be seen in Figs. 7, 8, this symmetry does not exist in our model, and in particular, we observe that $\theta_{1}\neq\theta_{2}$ after reflection. In fact, when reflected, axons remain much more parallel to the interface than light rays do.

Next we examine the case of an axon growing on a substrate with a graded stiffness of the form

where the parameter $m_{\chi}$ is the stiffness gradient.

Firstly, as a simple illustration of growth on a graded substrate, we initiate multiple axons on the boundary of a disk in a linearly increasing stiffness field, mimicking the experimental setup in Koser et al. (2016). The growth pattern obtained matches qualitatively that seen experimentally, see Fig. 9, where axons are gradually deflected towards stiffer regions (upper parts of the domain).

More generally, the local traction exerted by the growth cone exhibits a maximum on the line $\mathcal{D}$ of optimal stiffness $y=0$. We will show that $\mathcal{D}$ is a locally attracting set for our system in the sense that nearby trajectories are attracted asymptotically to $\mathcal{D}$. This property is a consequence of the nonmonotonic dependence of the traction force $T$ upon the stiffness $\chi$.

We consider an axon travelling in a neighbourhood of $\mathcal{D}$ of typical size $R$. Hence, the growth cone is initially located at the position $\boldsymbol{r}=rR\boldsymbol{e}_{y}$, with $r$ of order one, and travelling nearly parallel to the axis with $\boldsymbol{t}=\boldsymbol{e}_{x}+\theta\boldsymbol{e}_{y}+\dots$ with $\theta\ll 1$. We obtain the behaviour of this axon by expanding the traction (31) about $\chi^{*}$, while treating $R$ and $\theta$ as small parameters:

where $\mathbb{P}\left({\varphi_{0}}\right)\vcentcolon=\int_{-\varphi_{0}}^{\varphi_{0}}\boldsymbol{{e}_{\varphi}}\otimes\boldsymbol{{e}_{\varphi}}\otimes\boldsymbol{{e}_{\varphi}}\,\left({{\mathrm{d}{\varphi}}/{\varphi_{0}}}\right),$ and $\mathbf{P}\left({\varphi_{0}}\right)$ is defined in (33). Equation (56) extends (32) to higher order, which is necessary for our analysis, since the leading-order deflection of $O(R^{2})$ vanishes at $\chi=\chi^{*}$. Of particular interest is the transverse component $F_{n}\vcentcolon=\boldsymbol{F}\cdot\boldsymbol{n}$, given by

Since $F_{n}=0$ when $\theta=r=0$, there exist trajectories along $\mathcal{D}$ (however, with a small drift appearing at higher orders and controlled by $T^{\prime\prime\prime}(\chi^{*})$). For $r\gg\theta$, the first term in the parentheses dominates, and, since $T^{\prime\prime}\left({\chi^{*}}\right)<0$, the deflection component has opposite sign of $r$ and grows with $r$ in magnitude. Thus we interpret this term as an elastic-like restoring force which pulls the axon back to the optimal stiffness line $\mathcal{D}$. Conversely, very close to $\mathcal{D}$, i.e. for $r\ll\theta$, the deflection is dominated by the second term and takes now the opposite sign of $\theta$, which, again, will tend to reorient the axon towards the horizontal. We conclude that, close to the line of optimal stiffness, durotaxis has a restoring, stabilising effect which tends to maintain the axon near this line.

Note, however, that the longitudinal component $F_{t}\vcentcolon=\boldsymbol{F}\cdot\boldsymbol{t}\gg F_{n}$ dominates in the balance of forces, thus there may exist, in principle, conditions for which durotaxis will not suffice to keep the axon in the vicinity of $\mathcal{D}$. Furthermore, it cannot be directly concluded whether the stable trajectory is attracting (exponential stability), or whether axons will oscillate around it and repeatedly overshoot it.

To further explore the dynamics of the axon near $\mathcal{D}$, we perform numerical simulations shown in Fig. 10. Firstly, Fig. 10(a) illustrates that, in the vicinity of $\mathcal{D}$ and for small initial inclination, axons are attracted and converge towards this stable manifold. Second, Fig. 10(b) shows that there may exist an escape angle above which durotaxis is insufficient to keep the axon closer to $\mathcal{D}$ despite the negative durotaxis in that region. Clearly, if the axon crosses perpendicularly the line of optimal stiffness, it will not curve back. Whether nearby directions also escape in an infinite domain is neither easy to establish nor directly relevant to the problem. Indeed, in a much larger domain, it is possible that such trajectories curve back to the line of optimal stiffness. However, in practice such axons have migrated away from the region of optimal stiffness and are subject to other stimuli.

In conclusion these examples illustrate that, even for a simple uniform gradient, the nonmonotonic, nonlinear dependence of the locomotory force upon stiffness implies the existence of preferential regions where axons are expected to converge. Such a mechanism may be key to understanding neuronal focalisation towards regions of specific stiffness.