Curvature-sensitive kinesin binding can explain microtubule ring formation and reveals chaotic dynamics in a mathematical model

Microtubules are filamentous tubular protein polymers which are essential for a range of cellular behaviour, and are generally straight over micron length scales. However, in some gliding assays, where microtubules move over a carpet of molecular motors, individual microtubules can also form tight arcs or rings, even in the absence of crosslinking proteins. Understanding this phenomenon may provide important explanations for similar highly curved microtubules which can be found in nerve cells undergoing neurodegeneration. We propose a model for gliding assays where the kinesins moving the microtubules over the surface induce ring formation through differential binding, substantiated by recent findings that a mutant version of the motor protein kinesin applied in solution is able to lock-in microtubule curvature. For certain parameter regimes, our model predicts that both straight and curved microtubules can exist simultaneously as stable steady-states, as has been seen experimentally. Additionally, unsteady solutions are found, where a wave of differential binding propagates down the microtubule as it glides across the surface, which can lead to chaotic motion. Whilst this model explains two-dimensional microtubule behaviour in an experimental gliding assay, it has the potential to be adapted to explain pathological curling in nerve cells.


Introduction
The skeleton of cells (cytoskeleton) is essential for cell structure, dynamics and function.It is formed by long filamentous protein polymers of three different classes: actin, intermediate filaments and microtubules (MTs).Of these, MTs are the stiffest filaments with important roles in cellular processes, such as cell motility, division, organisation, adhesion, signalling and intracellular transport.
MTs are composed of α-and β-tubulin heterodimers which are bonded in a polar head-to-tail fashion to form long chains known as protofilaments; these protofilaments are then assembled into a helical tube.For a detailed description of how microtubules behave, see for example Hawkins et al. (2010) or Barsegov et al. (2017).
Many different proteins bind to MTs, controlling MT behaviours, including their nucleation, (de)-polymerisation, stabilisation, severing, biochemical modification, and crosslinking to each other or other cellular components (Lawson and Salas, 2013;Prokop, 2013).One particular class are MT-associated motor proteins, which use ATP as an energy source to walk along MTs, either to slide them against each other or to use MTs as intracellular highways to transport cargo around cells.Two fundamentally different classes of MT-associated motor proteins exist: the various members of the kinesin family of which most walk towards one end of the MT, and the dynein/dynactin complex which moves towards the other (Prokop, 2013;Schliwa and Woehlke, 2003).
Outside of cells, a powerful in vitro tool to study MT behaviour is a gliding (or motility) assay.In these experiments, motor proteins (typically kinesin-1) are adsorbed onto a solid surface in a drop of solution.When MTs are added, the surface-attached motor proteins attempt to walk along them, causing the MTs to glide over the surface.Typically in these assays, MTs stay relatively straight, as would be expected from their large persistence length (2 mm to 4 mm (Howard, 2001)).However, in certain experimental conditions, MTs can form micron-sized rings; such conditions include high MT density or the presence of an air-medium interface (Weiss et al., 1991;Amos and Amos, 1991;Liu et al., 2011;Kawamura et al., 2008;Kabir et al., 2012).Strikingly, these MTs are able to transform from straight gliding to a curved circling motion and back again (Liu et al., 2011), showing a dynamic and reversible ability to change curvature, implying that this is not due to permanent damage or irreversible damage/repair cycles (Schaedel et al., 2015) (see Figure 1).
Studying the mechanisms that underlie MT curling has important applications.For example, systems based on MT-kinesin gliding assays have potential uses as lab-on-a-chip medical devices, utilizing the ability to bind only selected proteins to MTs through the choice of specific cargo adapters, leading to advective transport rather than mere diffusion (Bachand et al., 2014;Chaudhuri et al., 2017).These nano-devices need to be robust for potential clinical uses, but the presence of MT rings may disrupt their design.Furthermore, curved MTs as observed in gliding assays are similarly found in cells, particularly in axons.Axons are the cable-like extensions of nerve cells; their structural backbone is formed by straight, parallel bundles of MTs.However, in the ageing brain or in nerves affected by certain neurodegenerative diseases (e.g.some forms of motor neuron disease), MTs are found to curl up with similar diameters as observed in gliding assays (Sanchez-Soriano et al., 2009;Voelzmann et al., 2016Voelzmann et al., , 2017)).
To explain this phenomenon, the model of local axon homeostasis has been put forward (Voelzmann et al., 2016).It proposes that MTs in axonal environments have a strong tendency to curl up likely due to high abundance of MTs and MT-associated motor proteins, thus meeting the conditions known to cause rings in gliding assays.Various MT-regulating proteins are required to 'tame' MTs into ordered bundles; functional loss of these regulators increases the risk of MT curling and could explain neurodegeneration linked to them (Voelzmann et al., 2016).This model represents a paradigm shift for the explanation of certain forms of axon degeneration, by putting the emphasis on MTs as the key drivers of axon decay.To lend credibility to this model, it is pivotal to identify and validate the mechanisms that can explain the phenomenon of MT curling.So far, Ziebert et al. (2015) introduced a model to explain the formation of MT rings which suggests that, in the presence of the MT-stabilising drug taxol, each tubulin dimer may exist in two distinct conformations, one slightly shorter than the other.In their model, protofilaments are able to switch between these two states; when only some of the protofilaments are switched this leads to a longitudinally curved MT as an energetically favoured condition, providing a mechanism to create rings via an internal change to the MT.
Here, we explore the complementary possibility that differential binding of external factors can actively contribute to MT curling.Peet et al. (2018) show that MTs which are being bent in a flow chamber normally straighten after the flow is removed, but stay curved in the presence of a nonmotile version of kinesin-1.They propose that this non-motile kinesin has a tendency to bind preferentially to the convex side of curved MTs and, by doing so, stabilise them in bent confirmation (see Figure 2); at higher concentrations this behaviour disappears, presumably because oversaturation occurs so the kinesin binds in equal amounts on all sides of the MT.
This behaviour is consistent with findings for other MT-associated proteins, in particular tau (Samsonov et al., 2004) and doublecortin (Bechstedt et al., 2014;Ettinger et al., 2016), which bind differentially between straight and curved MTs due to conformational changes that happen on the structural scale of the individual tubulin dimer: at a curvature of 1 µm −1 the tubulin dimer spacing at the outside of the MT is 2.5% larger compared to that of the inside.Here we present a model based on the hypothesis that curvature-selective binding can occur in MT-kinesin gliding assays; the flexible neck linker of the surface-attached kinesins can extend up to 45 nm from the surface (Palacci et al., 2016), and is therefore long enough to reach the curved sides of the MT which are typically held at around 17 nm above the surface (Kerssemakers et al., 2006)).Our model reproduces key behaviours of MTs observed in gliding assays, with a bistable regime where straight MTs and MT rings can coexist, and predicts how they can be controlled.Furthermore, we find unsteady propagating wave solutions and chaotic dynamics within the system, which have not been previously reported for filaments and may reflect true MT behaviours that have escaped the attention of experimenters so far.

Filament Dynamics
We consider a MT as an inextensible filament represented by a curve, x(s, t), parametrized by its arclength s measured at time t, which lies within a two-dimensional plane.This is a reasonable approximation for gliding assays, where the MTs remain close to the surface throughout their movement.The treatment shown here follows the standard case where the filament has no reference curvature (see, for example, Audoly, 2015;Ziebert et al., 2015;De Canio et al., 2017).
As the individual tubulin dimers are of fixed length, the total arclength is constant and we impose the inextensibility constraint, x • x = 1, where a prime denotes differentiation with respect to s. Relative to fixed Cartesian coordinate axes ex, ey, we define the tangent vector as where θ(s) measures the angle between t and ex.The normal vector n is then given by the relation t = κn, which defines the curvature κ = θ .
We will assume that the filament has a variable reference curvature (to be specified later), κ(s, t), and that the mechanical energy of the MT is a function of the squared deviation of the curvature from the reference curvature, Here λ is a Lagrange multiplier enforcing the inextensibility constraint and B is the bending (flexural) modulus.Taking the variational derivative of (2), utilizing δκ = n • δx , we find where δ represents the variation of a quantity and n = −κt.The elastic force density, f , acting on an element of the MT is therefore given by subject to x • x = 1.Here f has units of force per unit length, and may be considered as the circumferentially averaged surface stress (Lindner and Shelley, 2015).The MT is immersed within a viscous medium at very low Reynolds number, and so we use resistive force theory, a Stokes-flow approximation which takes advantage of the aspect ratio being small, = h/L 1, where h is the MT diameter (25 nm).This is the simplest approximation of slender-body theory, and gives the local dynamic relation (Lindner and Shelley, 2015), where v = ∂x ∂t = ẋ is the velocity of material points, f ext is the external force per unit length acting along the MT, µ is the fluid viscosity, and the tensor P ≡ (I + (ξ − 1)tt) = nn + ξtt reflects the anisotropic drag on the filament due to its shape.The constant c = ln(2 −1 ) is a free-space slender body ratio (Becker and Shelley, 2001), and we also use the free-space approximation ξ = 2 for an idealised slender filament.Both of these neglect the effect of the nearby surface; a more refined approach is likely to predict higher values of overall drag.There is nothing to prevent selfintersection of the filament within this model; if self-intersection does occur the filament is therefore assumed to go out of plane, crossing over or under itself.The equations of motion are therefore given by Due to the constraint, this is a ninth-order (in s) system of differential algebraic equations with index 3, so six spatial boundary conditions are required to fully specify the system.It is easier to work in terms of intrinsic coordinates which move with the filament, so we take the derivative of ( 6) with respect to s, and use ẋ = θn to give two equations in the normal and tangential directions respectively, where K = κ − κ is the excess curvature and τ = λ + BκK is a generalised tension.Here, as in Ziebert et al. (2015), we have assumed that the external force comes solely from the action of the kinesin motors, which force the microtubule along its tangent and so f ext = fmt, where fm is a constant.Equation ( 6) nevertheless allows the MT to move normal to its centreline.The position x may then be found from the filament angle θ by integrating (1).The natural boundary conditions for a free end of the MT come from the variational principle (3), and are given by For a fixed end we have x = x 0 , and therefore we set ẋ = 0 in (6) which leads to the force conditions, which are supplemented with either K = 0 for a freely rotating pinned end or θ = θ 0 for a clamped end.
For a straight filament with κ = κ = 0, (6) gives which allows us to estimate an appropriate tangential force being applied from measurements of the MT velocity.Typical velocities for gliding assays are 0.5 µm s −1 to 1 µm s −1 , with differences due to factors such as viscosity, ATP concentration, temperature and salt concentrations.While it would be reasonable to expect that the speed of the MT would be proportional to the amount of kinesin attached to the MT (and hence the surface kinesin density), as equation ( 10) implies, this is not seen in experiments, which show that the velocity is constant for kinesin densities of 10 µm −2 to 10 000 µm −2 (Howard et al., 1989).We therefore will use v to set an appropriate choice of fm.

Kinesin Binding
We now turn to the binding of the surface-bound kinesin to the MT, which we will model as a continuous field, assuming that the concentration is sufficiently high for this to be valid.Here we focus on the 'sides' of the filament, arbitrarily denoting them with + and − with associated bound concentrations c + and c − ; we do not model the protofilaments which are directly above the surface as we assume that they will not affect the reference curvature.The proteins bind and unbind to the filament according to standard protein binding kinetics, where a ± (κ) is a curvature-dependent association rate, δ is a disassociation rate, cmax is the maximum number of binding sites per unit length and D is a diffusion constant (measured as 0.036 µm 2 s −1 for kinesin-1 (Lu et al., 2009)).The left hand side of ( 11) is a material derivative, incorporating the fact that we are working in intrinsic coordinates while the kinesin is fixed to the surface, where vs is the instantaneous MT velocity (6) projected in the tangential direction, Dividing ( 11) by cmax, we use the bound ratios φ ± = c ± /cmax as dependent variables, giving At the ends of the MT, we allow no diffusion-based flux of the protein (although it will 'fall off' the trailing end with the velocity vs) and hence impose ∂φ ± ∂s = 0 at both ends.Although we are modelling the MT as a one-dimensional rod, in reality it has a complex protein structure, with each tubulin monomer consisting of approximately 450 amino acids folded into a 3D arrangement with charged residues protruding from the surface.Bending the entire filament moves these residues in relation to each other, expanding those on one side and contracting those on the other; such changes can be expected to change the binding kinetics of associated proteins, as these are also complex charged structures.The precise nature of this relationship is unknown, but here we will assume a sigmoidal relationship of the protein association rate a on the local curvature κ, where β acts as a scaling factor (with dimensions of length) to determine the degree of preferential binding of the kinesin to curved MTs.This implies that when the MT becomes curved the binding rates to each side of the MT will locally change, and increasing the value of β will mean that the differential binding is more sensitive to small curvatures.
The choice of this sigmoidal relationship ensures that the association rate both saturates at high curvature and that a ± (κ) is always positive; we have checked that other functional forms can be used to similar effect.
We assume that the average on-rate α 0 is proportional to both the number of kinesin molecules available in the vicinity of the MT and their ability to reach one side of the MT, where Γ is the surface kinesin density, assumed to be sufficiently large for depletion not to be a concern, d is the maximum distance kinesin can extend (45 nm (Palacci et al., 2016)), and ωon is an attachment rate per kinesin molecule within range, given as 20 s −1 (Chaudhuri and Chaudhuri, 2016).We note that this is a high estimate, because we neglect the binding to protofilaments directly above the surface.
Our final model assumption is that the local concentration of bound protein influences the intrinsic curvature of the filament, by acting as a brace on the side of the filament or some other conformational change, as suggested by Peet et al. (2018).As the protofilaments are bonded to each other via lateral bonds, it is assumed that this is able to affect the entire MT.If only one side of the MT has a high concentration this will prevent the filament from straightening, altering the MT reference curvature, while if both sides have bound protein then there will be no net effect on the curvature.Again, we assume a sigmoidal dependence of κ on the difference between the two 'sides' of the MT, where γ > 0 is a scaling factor that controls the steepness of the MT response to differential binding and κc > 0 is the maximum characteristic curvature.These unknown constants will depend on the precise nature of the bracing effect, but the measured lattice expansion of 1.6% in Peet et al. (2018) suggests κc = 0.625 µm −1 .The MT will therefore curve towards the side with less bound protein, as shown in Figure 2.

Non-dimensionalisation
We non-dimensionalise equations (7,12,13,16) with respect to the MT length L and the unbinding time δ −1 , resulting in the domain of integration being s ∈ (0, 1), yielding where the three dimensionless numbers, are the ratio of forcing to bending rigidity, called the flexure number in Isele-Holder et al. (2015), the ratio of elastic to viscous forces (inversely related to the Sperm number (Lowe, 2003)) and the updated base on-rate respectively.The boundary conditions at a free end are while at a pinned end we have, To connect back to the spatial positions, (6) may be written as, which we can calculate after solving for κ.We can also use (1) to get the shape, supplementing with (21) to find the position at a single point.
For the examples shown here, we will set D = 0.036 µm 2 s −1 (Lu et al., 2009), v = 0.5 µm s −1 , κc = 1 µm −1 (comparable to the 0.625 µm −1 suggested in Peet et al. (2018); the exact value does not affect the primary conclusions here), δ = 1 s −1 (Chaudhuri and Chaudhuri, 2016).For calculating χ, there is a wide range of measured values of the bending modulus B, depending on the measuring technique and conditions, with a noticeable length dependence (Pampaloni et al., 2006).Similarly, the viscosity µ is not clear, as the medium of the gliding assay is more viscous than pure water, and so we shall set χ = χ 0 /L 4 , where χ 0 = 1 µm 4 or 3 µm 4 in the examples below.
As we are considering only the protofilaments on the sides of the MT, we have cmax = 250 µm −1 if we only consider two protofilaments as being available for binding on each side.Combined with the estimates of the other values above, this gives α = 3.6 when Γ = 1000 µm −2 .

Uniform Steady-States
First we consider steady-state solutions with uniform-curvature, where κ = κ = κ 0 , τ , φ + and φ − are all constant.Evaluating ( 17) leads to values for the steady-state protein concentrations, and the following transcendental equation for κ 0 : The straight MT with no differential binding is always a solution to (23), while non-zero roots of κ 0 correspond to curved states with radius of curvature κ −1 0 .Note that (23) depends only on the parameters involved in the protein binding, not those connected to the mechanical response.Defining z = κ 0 /κc, b = βκc, steady-states are associated with the roots of the following threeparameter equation, defining the function j(z).Non-zero roots of j(z) will therefore correspond to steady-states with a non-zero uniform curvature, and so we wish to understand how the parameters affect the existence of these roots.If j (0) < 0, then a positive root must exist as lim z→±∞ j(z) = ±∞, which leads to a bifurcation condition from where non-zero roots intersect z = 0, As shown in Figure 3A, ( 25) defines two values of the binding rate α, (α 1 , α 2 ), at which pitchfork bifurcations occur; we focus on α as it is the parameter most available for experimental control.
A linear stability analysis (see Appendix) shows that the uniform solution becomes unstable with a single real positive eigenvalue at α = α 1 , until it re-stabilizes at α = α 2 .The first bifurcation at α = α 1 is supercritical, producing stable non-uniform solutions, but the solutions arising from α = α 2 can be either stable or unstable, depending on the values of b and γ, with a saddle-node bifurcation occurring further along the branch to re-stabilize solutions when they are initially unstable (Figure 3B).This saddle-node bifurcation occurs for a wide range of parameters and leads to multiple non-zero roots as shown in Figure 3A.These additional roots appear when j(z) is initially smaller than z but grows sufficiently fast to exceed z, giving a second non-zero root.
Provided that α 1 and α 2 are far enough apart, we also find that a nested series of Hopf bifurcations occur along both the unstable branches as the periodically spaced complex eigenvalues move across the real axis, generating unstable limit cycles (due to the positive eigenvalue still being present), with a corresponding reverse bifurcation when they pass back over the real axis.The positions of these Hopf bifurcations are indicated as points in Figure 3A and as curves on Figure 3B, which shows the parameter space as α and b are varied for a fixed γ = 1.A similar picture is found as the curvature-binding parameter γ is changed, but with a negative relationship between the two feedback parameters b (or β) and γ; when one is small the other needs to be large in order for the feedback strength to be large enough to create non-zero roots, as can be seen in ( 25).
We have therefore shown there exists a range of parameter values for which multiple stable steady-states exist, allowing for the simultaneous existence of straight and curved MTs at the same parameter values and for a single MT to be transferred between the two when suitably perturbed, as was shown experimentally by Liu et al. (2011) (Figure 1).As can be seen in Figure 3A, the non-zero constant value of κ 0 is generally below the prescribed characteristic curvature κc (i.e.z is less than 1); this value is approached for some parameter values but can be significantly less, allowing for variation in the ring sizes with a fixed κc.Furthermore, the Hopf bifurcations point to the potential existence of oscillating states, which we will explore below.

Numerical Solutions
We now solve the full set of partial differential equations, ( 17)-( 21), to see how the curvature of the MT evolves in time.To do this, we use the method of lines, discretising in s with a fourth-order central-difference formula to give a set of coupled nonlinear ordinary differential equations in t, which are solved using Mathematica.
As the straight solution is always a steady-state, we need to introduce an initial curvature perturbation to the MT to be able to see other behaviours.One option is to temporarily pin the leading tip of a gliding MT (as is observed in gliding assays when the MT encounters defective motors (Bourdieu et al., 1995)).For a large enough F, this will cause the MT to buckle and rotate around the pinned point (Sekimoto et al., 1995;Chelakkot et al., 2014;De Canio et al., 2017).We then allow the MT to unpin after some curvature (and therefore also differential binding) has been generated, and the MT will then either reach a curved configuration or re-straighten.This is shown in Figure 4 and Supplementary Movie S1, where two MTs that are unpinned after slightly different amounts of time settle into the two different steady-states.The resulting end-states correspond to the two stable steady-states seen at line III in Figure 3A.
Instead of pinning an initially straight MT, we can also generate an initial perturbation by bending the MT into a non-straight configuration and then allowing it to relax, mimicking an interaction with other MTs.In both of these cases, the exact basins of attraction of the two steady-states depends on both the binding and the mechanical parameters; a relatively large perturbation from straight gliding, which persists for long enough to produce binding differences, is required to move the MT into the curved configuration.

Propagating Waves and Chaotic Motion
As well as the steady-states detailed in Section 3.1, we also find parameter regimes where stable propagating waves move down the MT, shown in Figure 5 and Supplementary Movie S2; the resulting MT shapes look strikingly similar to those of cilia beating.These MTs show globally directed motion in the same way as straight MTs, gliding across the surface, but with the addition of a periodic oscillation feeding back from the tip.These sustained oscillations are caused by the reference curvature being translated along the MT as it glides.The curvature of the MT trailing tip, κ(L), shows the regular periodicity where the waveform reaches the tip (Figure 5B), and we denote the time between zeroes of κ(L) by T .
As the parameters are changed, the period and amplitude between these oscillations changes smoothly.Additionally, as the diffusion D is increased from D i = 0.036 µm 2 s −1 (Lu et al., 2009), we also see period-doubling bifurcations occur as shown in Figure 6, where T breaks symmetry, leading to a non-symmetric behaviour in the waveform (Figure 6C).As D is increased further a period-doubling cascade continues, leading to chaotic behaviour where the MT never settles into a periodic regime, as shown in Figure 6A.Upon further increasing of D, the wave solution disappears and the system settles into the uniform curvature solutions described in Section 3.1.While it is unlikely that D could be used as a control parameter in an experiment, these results demonstrate the range of possible outcomes in this system, and their sensitivity to parameter values.
After the initial period-doubling, the overall motion of the MT is biased, leading to the MT moving in a large circle whilst undulating, as shown in the insets of Figure 6A and in Supplementary Movies S3-S7.The radius of this large circle, Rc, decreases monotonically with D, as the waveform becomes increasingly asymmetric, shown in Figure 6D.
These oscillatory solutions, and the period doubling cascade to chaos, appear to exist in regions of parameter space around the second pitchfork bifurcation, α 2 , shown in Figure 3, provided that the 'wings' of the bifurcation diagram are large enough to include the Hopf bifurcations.They also exist when different sigmoidal functions κ = κcG(γ(φ + − φ − )) are used instead of ( 16), for instance where erf is the error function (results not shown).

Discussion
We have presented here a model for the occurrence of rings and arcs seen in MT gliding assays, based on the extrinsic effect of differential binding of the surface-attached kinesins which move the MTs forward.Our model goes beyond that of Ziebert et al. (2015), who assume that the formation of these rings is the result of an intrinsic property of the MTs.Our model incorporates the presence of the kinesins as external factors that drive a feedback loop where MT bending is recognised and stabilised, thus recruiting further kinesins; intrinsic properties of MTs as considered by Ziebert et al. (2015) may contribute to these effects and can be considered in our model by modifying the equation for the reference curvature ( 16).
Both models are able to reproduce the key aspect of the experiments, with regimes of bistability where both curved and straight MTs exist simultaneously for the same parameters, differentiated by the local forcing on the MT.In order for rings to be formed, we find that the system needs to be within a parameter regime (α, β, γ) which permits multiple steady-states, and the MT must undergo high induced curvature.Going beyond the model of Ziebert et al. (2015), our approach leads to two predictions which can be experimentally tested: 1.Our model predicts that varying the effective binding rate α should affect whether or not MT rings can form, as well as their size.This on-rate includes both the binding distance d, which can be experimentally altered by truncations of the kinesin-1 tail region as in the experiments performed in VanDelinder et al. (2016), as well as the surface-bound kinesin density Γ which may be directly varied by altering the kinesin concentration.2. Our model predicts that high MT density will encourage the formation of rings, consistent with the experiments of Liu et al. (2011).MT density contributes in two ways: it encourages the bending of MTs through MT-MT interactions, and it promotes the pushing down of MTs towards the surface of the assay during cross-over events, significantly enhancing the access to the convex MT side.
The non-kinesin parameters used in our model are already well-known or easily measurable, but the parameters β, γ and κc which connect the protein binding to the preferred curvature are entirely unknown.However, they could be characterised via fluid-bending experiments of the kind performed in Peet et al. (2018), and then fitting an appropriate variation of the model described here.The measurements of the curvature sensitivity for the protein doublecortin in taxol-stabilised MTs (Figure 3G of Ettinger et al. (2016)) suggest a value of β of approximately 3 µm.
Additionally, the exact values for the parameters involved in the combined kinesin on-rate α are not very well characterised.For instance, the surface kinesin density Γ is not routinely measured in gliding assays, but can be obtained via landing-rate experiments (Katira et al., 2007).Furthermore, other motor proteins than kinesin-1 (e.g.other members of the kinesin superfamily or the minusend directed dynein/dynactin motor complex), are able to translocate MTs in gliding assays.If parameter changes in these assays facilitate the occurrence of rings, this would provide additional data sets that could be compared and help to refine parameter determination.
The unsteady solutions where the MT oscillates while moving across the surface, including the period-doubling cascade to chaos, are particularly interesting mathematically as they are not immediately obvious from the governing equations.An oscillating regime for a clamped or pinned filament (without the preferred curvature) driven by the tangential motors, as considered here, was found by Sekimoto et al. (1995), with a Hopf bifurcation occurring above a critical forcing F; De Canio et al. ( 2017) also show this behaviour for the case of a filament with a follower force acting on the free end, and expand upon its origin.The behaviour seen here is similar, with a self-sustained oscillatory waveform, but these authors did not report chaotic dynamics.These oscillations may be biologically relevant; the wavy MTs generated by our model look similar to those shown in an experimental figure of Gosselin et al. (2016), as well as a MT seen in Supplementary Movie 1 of Scharrel et al. (2014), and similar regularly undulating curved MTs are seen in cells (Brangwynne et al., 2007(Brangwynne et al., , 2006)); although these are assumed to be caused by mechanical buckling, it is possible that this kind of curvature-dependent binding may enhance the effect.There may also be a connection to the MT phenomenon described as 'fishtailing', where MTs oscillate laterally whilst their head is stuck, as shown in Applewhite et al. (2010) and Weiss et al. (1991) for example.We therefore encourage experimentalists to look out for this kind of behaviour.

Future Directions
Our current model is only a starting point, and a number of further aspects can or should be incorporated.
1.The arrangement of protofilaments into the MT is via a helical arrangement, with a skew angle inducing a global supertwist for the MTs where moving forward along one protofilament involves rotating around the MT; the exception is for MTs with 13 protofilaments which are almost straight, and this is therefore the type considered in our model so far.MTs with more or fewer protofilaments are shown to rotate as the kinesin moves along them (Ray et al., 1993), inducing a torque that could be considered in future versions of the model, particularly as axonal MTs are not always the 13 protofilament type (Chaaban and Brouhard, 2017).Furthermore, Kawamura et al. (2008) find that more rings occur in gliding assays when they use freshly prepared MTs, which have fewer 13-protofilament MTs than their 24 hour aged MTs, suggesting that this MT rotation might enhance the ring-formation by making it easier for the kinesin to be bound to the outer side of the MTs; this effect may be explained by kinesin stepping from one protofilament to another as it reaches its maximum extension, as it does to move around obstacles (Schneider et al., 2015).
2. The inclusion of the crosslinking protein streptavidin into gliding assays induces the formation of bundles of curved MTs, known as spools, where multiple MTs are attached together and the entire structure is bent into an arc (VanDelinder et al., 2016).Luria et al. (2011) found more smallcircumference spools than predicted by their simulations; these spools have a similar diameter to that of the non-crosslinked MT loops (Liu et al., 2011;Kawamura et al., 2008), and we suggest that the mechanism we propose may be responsible.In particular, Lam et al. (2014) found that the size and number of crosslinked spools depends upon the kinesin density, with the presence of more kinesin leading to fewer small-diameter spools, which is consistent with our results.
Incorporation of MT-MT interactions via crosslinking proteins would therefore be a natural extension of our model.
3. The core of this model is suitable for other situations where differential protein binding may influence filament curvature, for instance where the protein is freely diffusing in solution as in Peet et al. (2018).The model assumes that the MT stays in the same vertical plane, as it is attached to the surface by the kinesin, and the extension to three dimensions may be required to properly model the situation in cells, incorporating both the twisting as well as the bending of the MT, as well as modelling all the protofilaments individually.
4. As mentioned in the introduction, other MT-binding proteins occurring in nerve cells, such as tau and doublecortin, have been shown to be curvature-sensitive, and we therefore recommend that these are tested to see if they can also reinforce the curvature.Additionally, an extension to the model to incorporate competition between proteins for the same binding sites may explain how certain proteins have a protective function.

Fig. 1
Fig. 1 Time-series of an initially straight MT forming a loop at 35 s, rotating until 75 s before re-straightening, extracted from Supplementary Movie 1 of Liu et al. (2011).Another adjacent MT stays in a loop through the whole video.Each frame is 15 µm square.Used with permission of Prof. J. Ross.

Fig. 2
Fig.2Sketch of the proposed model.A) When the MT is straight, the surface-bound kinesin is equally likely to bind to each side, with no effect on the overall curvature.B) If a MT becomes curved, the likelihood of kinesin binding from each side becomes asymmetric; this asymmetry in the amount of bound kinesin to each side induces curvature by acting as a lateral reinforcement.As the MT then continues to glide across the array, the newly encountered kinesin will also preferentially bind to the same side.

Fig. 3
Fig. 3 (A) Bifurcation diagram for b = 3, γ = 1, plotting normalised curvature z against dimensionless binding rate α.Unstable branches are marked by dashes, dots represent the location of the Hopf bifurcations shown in B. Vertical lines I-IV show positions with 1, 3, 5 and 1 values of κ 0 respectively.(B) Parameter space map (over binding rate α and preferential binding affinity b) showing the position of the three types of local bifurcations seen in the system, for γ = 1, χ 0 = 3 µm 4 , L = 10 µm.The dotted line corresponds to the bifurcation diagram shown in (A).

Fig. 5
Fig. 5 A) Shape evolution of a MT showing oscillatory behaviour as it moves across the surface.At these parameter values propagating waves in κ cause the MT to periodically undulate as it moves across the surface, settling into a steady rhythm.One period of these shapes are shown superimposed in (B).For these parameters, the trailing tip preferred curvature, κ(L, t), plotted in (C), has a regular gap between zeroes of T = 3.466.Here α = 8, b = 3, γ = 1, χ 0 = 1 µm 4 , L = 5 µm.

Fig. 6
Fig. 6 (A) Time T between the zeroes of the trailing tip curvature κ(L, t) (after initial transient behaviour) as the diffusion D is increased, showing a period doubling cascade to chaos.Inset above are plots (all on the same scale) showing the motion of the trailing tip of the MT after the initial transient behaviour.For D/D i = 3.55, 3.56, 3.57, the MT settles into a steady orbit with decreasing radius, whereas for D/D i = 3.5823 the MT moves around the plane in a chaotic manner.(B) Superimposed shapes for D/D i = 3.57 over a sub-interval of the orbit, showing a larger amplitude than in Figure 5. (C) Plot of the reference curvature at the end of the MT, κ(L), showing the two values of T which can be seen in (A).(D) Plot showing how the radius of the large orbit traced out by the oscillating MT decreases monotonically as D/D i is increased.All other parameters are the same as in Figure 5.